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1.  PROBLEM  DESCRIPTION 

V 

Many  important  problems  in  fluid  dynamics,  among  other  areas, 
are  modeled  by  nonlinear  parabolic  differential  systems  with  initial 
values  given  in  one  ^independent  variable*  x,  and  boundary  values  in 
the  remaining  dependent  variables.  Hyperbolic  systems  can  sometimes 
be  treated  as  a  special  case.  For  example,  the  inviscid  flow  case 
of  the  Navier-Stokes  equations  is  a  hyperbolic  system,  which  the 
viscous  flow  case  is  elliptical.  A  survey  of  currently  used  numerical 

methods  is  in  Richtmeyer  and  Morton  [2].  In  subsonic  flow  cases,  the 

✓ 

nonlinear  terms  are  small  enough  to  be  Ignored,  but  these  terms  must  be 
included  in  supersonic  and  hypersonic  flow.  These  numerical  calculations 
usually  Involve  a  finite  difference  mesh  over  the  boundary  value  problem 
variables,  resulting  in  a  space  discretization  matrix  equation  which  for 
the  nonlinear  system  varies  at  each  step  in  x,  the  independent  variable 
representing  time  in  the  dynamic  case  or  one  of  the  space  variable  for 
the  steady  state  case.  Then  this  nonlinear  system  is  solved  as  an 
initial  value  problem  in  x.  The  initial  value  problem  is  usually  solved 
by  a  one  step  implicit  method  for  reasons  of  cost  and  stability.  Some 
methods  based  on  finite  element  methods  for  the  boundary  value  problem 
can  be  used,  but  successful  methods  are  only  available  for  the  linear 
cases,  such  as  subsonic  flow  problems  - 

All  of  these  methods  require  large  amounts  of  computer  memory  to 
store  the  matrices,  and,  particularly  in  the  nonlinear  case  where  the 
matrices  must  be  reevaluated  often,  large  amounts  of  time.  Therefore,  it 
is  desirable  to  investigate  the  relationship  of  various  aspects  of  these 


->Wfi  Vim'  i*i *■*£&*?■■&*  '  •  •*'*«&* 


numerical  methods  in  an  effort  to  reduce  the  total  computation  time  with 
no  loss  in  accuracy  or  significant  increase  in  storage  requirements.  To 
give  an  overview  of  the  current  state  of  development,  a  sample  problem 
which  has  been  studied  by  the  principal  investigator  [4]  is  described. 

The  incompressible  fluid  flow  around  a  cone  at  hypersonic  speed 
and  angle  of  attack  a  >  0  is  modeled  by  a  parabolic  system  of  nonlinear 
partial  differential  equations  expressing  conservation  of  energy,  mass, 
and  momentum,  plus  an  algebraic  equation  of  state.  Typical  flow  variables 
are  functions  of  density,  velocity  and  energy.  The  asymptotic  (steady- 
state)  solution  in  three  dimensions  is  sought.  A  suitable  coordinate 
system  for  a  cone  shaped  object  uses  the  variables  x,  the  length  from 
the  tip  along  the  cone  generator;  n  the  normal  to  the  surface  relative 
to  the  bow  shock  stand-off  distance  (n  =  £/d(x,<f)  where  £  is  perpendicular 
to  x  and  d  is  the  bow  shock  stand-off  distance  computed  from  theory) ;  and 
<|>  =  180°  at  the  leeward  side.  Separation  is  likely  to  occur  at  the 
leeward  side  at  significant  angle  of  attack  a  >  0,  and  standard  numerical 
methods  have  proven  inadequate  to  model  this  case,  so  special  computer 
methods  have  been  developed  for  it.  See  Figure  1. 

Lubard  and  Helliwell  [5]  have  treated  this  problem  as  a  parabolic 
boundary  value  problem  in  if  and  n  since  theoretical  results  are  available 
on  the  behavior  of  the  bow  shock,  and  as  an  initial  value  problem  in  x 
which  allows  a  marching  type  numerical  solution  to  be  generated  given 
an  initial  condition  away  from  the  point  x  ■  0.  They  treated  the  non¬ 
linear  system 

3U  3F(U)  3G(U)  3V(U,  3U/3n,  3U/3<f)  3W(U,  3U/3n.  3U/3<f) 

3x  3n  a*  =  3n  3<f  (1) 

where  U  is  the  m-dimensional  state  vector,  F  and  G  are  given  vector 
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functions  of  U,  and  V  and  W  are  vector  functions  of  U  and  its  partial 
derivatives.  Steady  state  Navier-Stokes  equations  are  a  special  case 
of  (1).  For  the  purpose  of  illustration,  let  the  right  hand  side  of  Cl) 
be  0,  giving  the  equation  for  inviscid  flow.  A  discussion  of  boundary 
and  initial  conditions  can  be  found  in  [5].  By  using  a  finite  difference 
scheme  in  n  and  <j>  for  the  boundary  value  problem,  a  one-step  implicit 
integration  scheme  in  x  can  be  employed  to  solve  the  initial  value 
problem.  For  the  implementation  of  Lubard  and  Helliwell  [5],  central 
differences  at  each  x  value 


S3.  =  JL/i, 

0 An  f  U 


Hi  =  2Ah  (U(Vl*  V  *  U(Vi’  Vy 


’If  '  ^  (B(V  *>-l>  -  2U<V  +  0<V  Vi*) 


(2) 

(3) 


are  used  with  the  Backward  Euler  implicit  formula 

U^xi+1’  nj*  ^  =  U^xi’  nj  ’  +  Ax  8U(xl+l’  nj‘  (t,k)/8x*  (*) 

Stability  considerations  derived  from  considering  the  numerical  solution 
of  an  associated  linearized  system  of  equations  leads  to  both  lower  and 
upper  bounds  on  Ax  as  a  function  of  An  and  A$. 

To  understand  the  implementation,  consider  the  linearized  problem 
for  A  =  3F/3U  and  B  *  3G/3U  given  as 


,./ 3u(!W  5U(V  ! 

Applying  the  trapezoidal  rule  u(xi+1)  =  U(x^)  +  — ji  — ^ ^ — (  , 
a  more  accurate  implicit  scheme  than  Backward  Euler,  and  using  a  truncated 


Taylor  Series  for  F(U)  and  G(U)  given  by 
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F(U(xi+l))  "  F±+1  “  f1  +  A1  (U(xi+1)  -  U(xt»  +  0(Ax2) 
G(U(xi+l))  =  °1+1  =  ci  +  B1  <U(X1+1>  ~  U<xi^  +  °^x2> 


yields  the  system  of  linear  equations 


I+f* 


3A 

3h 


3B 


i'\ 


3<t>  J 


..i+1 


+  *2.  f  1*1  +  10  U 

+  2  \3n  +  3* /Jj 


3G  ■ 

3*j 


This  large,  sparse  system  can  be  solved  by  methods  such  as  the  alternating 

direction  implicit  (ADI)  method  of  Douglas  [6]  which  solves  only  the 

equation  in  n  first,  then  the  equations  in  <j>.  Beam  and  Warming  [7] 

3 

introduce  an  error  of  (Ax)  in  an  approximate  factorization  scheme  based 
on  Peaceman  and  Rachford  [8]  by  replacing  (6)  with 


! 


l+~ 


Ax  3A 


2  3n 


\  (  i 

: !  I  +  *2. 

l!  2  3<f> 


-Ax 


3F 

an 


,i+l 


-  1  +  -F: 


Ax  3A' 


3B 

3<J> 


+  d- 


3gV  OZ-A  \ 

a *y  +  0(&x  >■ 


(7) 
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Since  the  error  introduced  is  of  the  same  order  Ax  as  the  error  in  the 
trapezoidal  rule,  stability  is  not  affected.  This  equation  can  then  be 
solved  in  two  levels 


This  method  has  several  disadvantages.  If  used  with  a  more 


accurate  difference  method,  the  error  introduced  in  the  factorization 
will  lower  the  error  order  of  the  method;  however,  if  used  with  a  lower 
order  method  such  as  the  Backward  Euler,  good  results  could  be  expected. 


However,  Lubard  and  Helliwell  note  that  the  difference  in  *  near  n  =  0 

has  a  singularity  there,  and  use  instead  a  method  that  uses  the  factorization 

(7)  in  n  but  solves  for  each  set  of  solutions  at  each  ^  in  sequence 

<b  =  0°  to  6  =  180°  in  steps  of  A<f>.  This  is  done  iteratively  until  the 
Y0  k 

computation  converges.  At  each  the  resulting  system  of  linear 

equations  is  an  n  by  n  block  tridiagonal  matrix  of  block  size  m  by  m 

where  i  =  1,  ...»  n  for  and  m  =  6,  the  number  of  states  at  each  point. 

The  method  is  comparable  to  a  Gauss-Seidel  iteration  with  each  element 

2 

of  the  solution  replaced  by  a  (n*m)  size  linear  equation. 

In  actual  practice,  the  equations  are  rewritten  to  compute  the 
change  AU  in  the  current  value  of  U(Xj+^).  This  is  called  the  delta 
form  of  the  corrector  and  yields  a  linear  block  tridiagonal  system 


B1  C1 


A2  B2  ^2 

*  .A  B 
n  n 


I _ 

where  A^,  B^,  and 
for  u(xi+i) • 


i 

I 

j  AU  =  RHS 


are  square  ra  by  m  matrices. 


(8) 


RHS  is  a  m  by  n  corrector 
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2.  Research  Conducted 

A  portable  program  called  HVSL  [9],  which  is  a  modified  version 
of  the  Lubard  and  Helliwell  code  and  is  used  at  Arnold  Engineering 
Development  Center  (AEDC),  was  available  for  experimentation.  It  is 
portable  since  the  initial  state  of  the  system  is  read  in,  in  part, 
from  cards  and  then,  after  solving  a  boundary  value  problem,  all  initial 
values  at  50  n  points  and  19  <j>  points  are  known.  These  can  be  changed 
by  an  interpolating  procedure  included  in  the  code.  For  experimental 
purposes,  the  number  of  tj>  values  from  windward  (0°)  to  leeward  (180°) 
was  reduced  from  19  to  3.  Validation  tests  were  run  at  angle  of  attack 
a  =  1°  to  determine  if  this  modified  system  produced  the  same  solution. 
At  least  three  decimal  place  agreement  was  observed  at  $  =  0°,  90°,  and 
180°,  so  it  was  concluded  that  this  much  less  expensive  test  program 
was  adequate  for  testing  modifications  to  HVSL. 


The  following  changes  were  made  to  HVSL. 

1.  The  initial  (predicted)  value  of  U(x^+1),  representing  the 
solution  U(x^+^,  nk>  for  k  =  1,  NK,  A  =  1,  NL  by  the  Euler  explicit 
method,  is  computed  by 


=  u(xj)  +  (U(xJ  “  u(xj_p)  *  Axj+i/Axj' 


While  this  appears  to  be  using  a  finite  difference  to  approximate  the 
derivative  U(x^)  =  dU(x^)/dx,  it  is  actually  the  correct  value  since,  if 
Euler's  implicit  formula  is  iterated  to  convergence, 


U(Xj)  =  U(x^_1)  +  AXj  U(Xj), 


The  modified  program  stores  the  derivative  term  Ax^U(Xj)  in  AU  after 
the  last  evaluation  and  correction  of  whatever  numerical  method  is  in  use, 


and  this  is  used  in  an  Euler  explicit  predictor.  No  additional  storage 


I 

A  I 


is  required. 

2.  The  check  for  convergence  in  subroutine  IMPETA  checks  to  see 
if  the  right  hand  side  of  the  matrix  equation  (8)  satisfies 


u  O  C. 

Z  RHS,  .  ,  <  6*10 

i=1  i.l.k.l  - 


for  all  q.  ,  for  each  <p.,  Z  =  1,  NL,  for  convergence  at  x  .  The  subscript  i 

K.  X/  J 

refers  to  the  six  state  variables.  In  test  runs,  no  calculation  ever 
terminated  due  to  meeting  this  convergence  test,  but  instead  the  maximum 
number  of  iterations  were  used.  A  more  appropriate  convergence  criterion 
would  be  to  stop  when  the  last  corrector  did  not  make  any  changes  in  the 
third  decimal  place  of  any  variable,  and  this  relative  change  criterion 
was  implemented. 

3.  The  Lubard  and  Helliwell  code  uses  the  Backward  Euler  corrector 


U(Xj+i>  =  U (Xj )  +  Ax  U(x^+1) 


to  calculate  successively  better  approximations  to  U(x^+^).  Using  the 
the  delta  form,  it  is  seen  that 


AlP  =  U° (x. , , )  -  U(x.)  =  Ax  U(x . ) 


m1  -  u1^)  -  u°(xJ+1) 

=  U(x.)  +  Ax  U^x.j,)  -  (U(x.)  +  Ax  U(x,)) 
J  3+1  J  j 


=  Ax(U (Xj+1)  -  U(Xj)) 
iU1  -  »*(xj+1)  - 

=  Ax(Ui_1(xj+1)  -  Ui_2(xj+1)),  i  -  2,  ....  5. 


I  " 


Since  this  calculation  is  already  programmed,  AU  can  be  used  as  is  in 


two  different  corrector  formulas.  The  Trapzoidal  Corrector 

U1(x)  =  U(x  )  +  Ax(Ui-1(x  )  +  U(x,))/2. 

J+l  j  j  3 

can  be  implemented  by 

ul(Xj+l>  =  U^xj^  +  ^^(x^)  +  U(x^))/2. 

=  U(Xj)  +  Ax  U(x^)  +  Ax(U^(Xj+^)  -  U(x^))/2. 

=  U°(xj+1)  +  AU1/2.  (13) 

and  then 

U1(xj+i)  =  U(x.)  +  Ax(Ui-1(xj+1)  +  U(Xj))/2. 

=  U(Xj)  +  Ax(Ui_2(xj+1)  +  U(Xj))/2  +  Ax(Ui_1(xj+1)  -  Ui_2(xj+1)/2 
=  U1_1(xj+1)  +  AU1/2.  (14) 

The  Iterated  Multistep  Method  (IMS)  due  to  Hyman  [10]  is: 

U^(Xj+^)  =  U(x^)  + 

yl^Xj+l^  =  U^Xj^  +  ^(^(x^)  +  U(x  ))/2. 

U1^  x)  =  U1_1(xj+1)  +  Ax(Ui_1(xj+1)  -  Ui-2(xj+1))/(i+l),  i  =  2,3,4,.... 

This  can  be  formulated  the  same  as  the  Trapezoidal  method  for  i  =  0,1, 
and  then 

Ui(xj+1)  =  U1-1(xj+1)  +  (AUi)/(i+l),  i  =  2,3,....  (15) 
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These  alternative  methods  have  both  stability  and  accuracy  advantages 
over  the  implicit  Euler  corrector.  The  Trapezoidal  corrector,  when 
applied  to  the  linear  complex  equation 


U  =  X  U,  (16) 

with  nonzero  initial  value  of  Uq,  damps  out  any  error  introduced  by 
either  machine  roundoff  error  or  the  discretization  of  the  solution  with 
respect  to  x  for  any  Ax  >  0  as  long  as  X  has  a  negative  real  part.  This 
is  called  A-stability  [11].  The  exact  solution  to  (16),  exp(Xx)Ug,  behaves 
the  same  way  since  an  initial  error  d^  yields  the  solution  exp(Xx)UQ  +  exp(Xx)dg, 
and  thus  the  error  contribution  goes  to  zero  as  x  goes  to  infinity  if 
X  has  negative  real  part.  This  assumes  that  the  Trapezoidal  corrector  is 
solved  exactly,  which  is  possible  in  the  linear  case  since 

U±+1  =  (1  +  AxX/2)  /  (1  -  AxX/2)  U±. 

Stability  behavior  is  somewhat  different  if  is  solved  iteratively, 

as  must  be  the  case  in  a  nonlinear  equation.  Thus,  the  stability  behavior 

of  the  trapezoidal  corrector  should  be  investigated  further. 

The  accuracy  of  the  Trapezoidal  corrector  is  based  on  the  local 

discretization  error,  which  is  the  size  of  the  error  in  if  were 

2 

the  correct  solution.  This  is  proportional  to  Ax  for  the  implicit 

3 

Euler  corrector,  but  to  Ax  for  the  Trapezoidal  corrector.  Thus,  the 
Trapezoidal  formula  is  more  accurate. 

The  IMS  method,  applied  to  (16),  has  the  property  that  each  successive 
corrector  iteration  increases  the  stability  region,  i.e.  AxX  such  that 
errors  introduced  are  not  increasing  in  size  as  the  calculation  proceeds. 

10 


Figure  2 


Linear  Stability  analysis  of  Iterated  multistep 
method.  Consecutively  larger  figures  are 

U1  for  i  -  0,1, 2, 3. 
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See  Figure  2.  Also,  in  the  linear  case  only,  each  application  of  the  IMS 
corrector  equation  increases  the  accuracy,  i.e.  the  discretization  error 
of  U*  is  proportional  to  Ax*+^.  In  the  nonlinear  case,  the  error  term  is 
similar  to  that  of  the  Trapezoidal  corrector. 

The  above  changes  were  made  to  the  HVSL  test  program,  and  the  resulting 
values  were  compared  to  the  original  program.  It  was  noted  that  none 
'j£  the  test  cases  achieved  convergence,  either  the  old  or  new  convergence 
criterion.  However,  all  significant  numerical  values  did  agree  to  two 
decimal  places,  so  it  was  concluded  that  there  was  a  marginal  stability 
problem,  and  the  stability  analysis  in  [12]  was  insufficient  to  explain 
the  phenomenon  since  that  analysis  was  based  on  linearizing  the  system 
and  inspecting  the  eigenvalues  of  the  resulting  Jacobian  matrices.  There¬ 
fore,  a  simpler  test  case  involving  only  one  space  variable  and  time 
as  the  independent  variable  was  used  to  study  the  three  methods.  The 
nonlinear  problem  has  properties  similar  to  the  HVSL  problem,  and  uses 
the  same  discretization  as  the  Lubard  and  Helliwell  method. 

The  quasi-one-dimensional  time  dependent  flow  of  an  inviscid  perfect 
gas  through  a  converging-diverging  nozzle  use  the  variables: 
x  =  distance,  normalized  to  [0.,!.], 

A(x)  =  nozzle  cross-sectional  area, 


H 


where  the  state  variables  are  P,  the  gas  density,  in  =  pu  where  U  is  the 
velocity  along  the  x  axis,  and  e  =  p(e  +  U2/2)  for  e  =  c  T,  where  T  is  the 

V 

temperature  and  cv  is  the  gas  constant.  Thus,  T  =  (e  -  m2/(2p))/(pcv> . 

The  equations  are 


3U  3F 
3t  3x 


G 


(17) 
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where 


F 


Ul  -  , 

pRT  +  m/p 

I 

|  n  e/p  +  mRT  j 


and 


G 


-  m  d(ln  A)/dx  ! 

-  m^  d(ln  A)/dx/p  j 

-  m(e/p  +  RT)  d(lnA)/dx  , 


A  working  version  of  this  program  was  provided  by  Dr.  John  C.  Adams 
of  AEDC.  The  program  linearized  (17)  with  respect  to  t,  giving 


dU  3F  dU 
dt  3U  dx 


(18) 


and  replacing  —  by  finite  difference  approximations  over  a  net  of  101  equally 
spaced  points.  Thus,  the  result  is  a  system  of  303  ordinary  differential 
equations  in  t,  with  algebraic  boundary  condition  consistent  with  the 
method  >f  characteristics  solution  in  one  dimension  at  the  inlet,  and 
extrapolation  of  supersonic  outflow  at  the  exit.  The  resulting  block  tri¬ 
diagonal  system  is 


A 


2 


AU  =  RHS 


B 

n 


for  3  by  3  matrices  A±,  B  ,  1=1,  ....  101.  The  linearized  initial 

value  problem  is  solved  exactly,  once  each  time  step.  The  numerical  method 
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is  parametrized  as  [13] 


d0<tn> 

dc 


1  (1  +  z)  A  -  zV 

At  1  +  0A 


0(tn) 


where  A  is  the  forward  difference  operator,  V  is  the  backward  difference 
operator,  and  0=1,  z  =  0  yields  the  exact  implicit  Euler  solution 
of  the  linearized  problem;  0=1/2,  z  =  0  yields  the  exact  trapezoidal 
solution;  and  0=0,  z  =  0  would  yield  the  explicit  Euler  predictor 
except  the  program  would  divide  by  zero  if  0  =  0.  This  does  not  emulate 
the  iterative  solution  technique  in  the  Lubard  and  Helliwell  code,  so 
the  program  was  rewritten  to  use  an  explicit  Euler  predictor  approximated 
by  0  =  10  r,  r  large,  then  successive  linearizations  and  computations  of 
AU*  consistent  with  the  implicit  Euler  technique  used  in  HVSL.  Then 
either  AU*  could  be  used  as  is  to  get  the  implicit  Euler  predictor,  or 
else  equations  (10,  13,  14)  for  Trapezoidal  corrector  or  (10,  13,  15)  for 
IMS  corrector  could  be  used. 

The  results  for  the  Trapezoidal  and  IMS  test  runs,  using  a  constant 
4  corrector  iterations ,  agreed  to  4  decimal  places  with  the  original 
program  at  the  4th  time  step,  and  to  2  decimal  places  after  15  time  steps 
(the  equations  are  being  integrated  to  steady  state).  The  velocity  U, 
which  depends  on  e  and  p,  becomes  unstable  by  the  4th  time  step  when  the 
implicit  Euler  corrector  was  iterated  4  times,  by  the  10th  time  step  when 
iterated  only  once.  Thus,  the  stability  properties  of  the  linearized 
equation  are  seen  to  be  different  from  those  of  the  nonlinear  equation.  Note 
that  the  step  size  At  was  chosen  to  meet  the  Courant-Friedrich-Lewy 
criteria  for  the  linearized  implicit  Euler  formula,  yet  this  step  size 
does  not  work  with  the  nonlinear  equation  it  is  based  on.  This  confirms 
that  a  stability  problem  exists. 
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A  new  technique  has  been  developed  to  study  stability  of  ordinary 
differential  equation  integrators  as  they  are  applied  to  nonlinear  differential 
systems  [14].  The  analysis  can  usually  be  carried  out  using  an  interactive 
graphics  package  called  STAN.  If  the  user  knows  an  approximate  equilibrium 
point  U*  where  dU/dt  ■  0,  then  it  is  possible  to  investigate  the  stability 
of  any  two  dimensional  subsystem  by  varying  only  two  values,  i.e., 
we  investigate  the  system  u  =  (u^,  u^)  =  U*  +  e^d^  +  e^d^  where  e^,  e^ 
are  the  i-th  and  j-th  unit  vectors  and  d^,  d^  are  scalar  perturbations. 

This  allows  the  contractivity  region  defined  by  Dahlquist  [15]  to  be 
mapped.  The  boundary  of  this  region  consists  of  points  at  which  the 
forward  difference  in  the  independent  variable  t  of  a  quadratic  function 
V(u)  =  u*Qu  is  zero.  V(u)  is  chosen  in  the  same  way  as  the  Liapunov 
stability  function  is  chosen  [16],  using  the  Jacobian  matrix  of  the 
derivative  with  respect  to  the  vector  u  =  (u^,  u^).  Thus,  if  this  boundary 
is  found  using  the  exact  solution,  then  any  solution  U(tg  +  At)  with 
initial  conditions  U(tQ)  on  the  boundary  has  the  property  V(u(tQ  +  At))  = 
V(u(tQ)).  This  can  be  further  refined  to  compute  a  stability  region 
inside  which  u(t^  +  nAt)  will  stay  for  any  n,  at  least  in  the  autonomous 
case.  Similar  regions  can  be  generated  for  the  numerical  solutions  using 
the  same  At.  Figure  (3)  illustrates  the  results  for  the  linear  case: 

y  =  u1  +  iu2 

y  =  Xy,  y (0)  =  (u^  u2)  (19) 

X  =  ln(yQ)/At 

with  exact  solution  at  mesh  points  y(t^)  =  yg+^.  F°r  v(u)  “  ui  +  u2’  t*ie 
contractivity  region  and  stability  region  are  both  the  unit  circle  about  0. 
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Figure  3A.  Non-Linear  Stability  analysis  of  y'  •  Xy: 

Contractivity  region  for  a)  explicit  Euler, 

b)  trapezoidal  with  one  corrector  iteration, 

c)  trapezoidal  with  2  corrector  iterations. 
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Non-Linear  Stability  Analysis  of  y  -Ay 
a)  explicit  Euler  predictor,  one  implicit 
Euler  corrector  b)  same  predictor,  2 
corrector  iterations. 


The  point  (0,0)  corresponds  to  the  point  at  infinity  in  the  linear  analysis 
described  earlier,  and  the  point  (1,0)  corresponds  to  the  origin  in  the 
linear  analysis. 

To  get  the  exact  solution  to  a  differential  system,  a  program 
SERIES  [17],  written  in  PASCAL,  is  used.  It  can  accept  up  to  20  first 
order  ordinary  differential  equations  as  input  and  will  write  a  FORTRAN 
subroutine  S0L(T,  Y,  YO,  IND)  which,  when  called  for  a  certain  value  of 
t  =  T,  given  initial  conditions  of  u(t)  at  t  =  0  stored  in  the  FORTRAN 
array  YO,  will  recursively  generate  the  coefficients  of  the  power  series 
solution  of  u(t)  starting  with  the  constant  terms  stored  in  YO,  and 
output  in  the  array  Y  the  solution  u(t) ,  or  else  indicate  that  the  radius 
of  conveigence  of  the  series  is  not  greater  than  t  by  setting  IND  to 
certain  nonzero  values.  Of  course,  not  every  possible  function  is 
included,  but  the  trigonometric,  logorithmic,  and  exponential  functions 
of  the  dependent  variables  are  allowed,  as  well  as  most  FORTRAN  functions 
of  the  independent  variable.  This  program  has  been  tested  on  numerous 
nonlinear  systems  and  the  resulting  subroutine  SOL  has  been  interfaced 
to  STAN,  but  no  system  of  partial  differential  equations  has  been  included. 

Two  different  stability  analyses  were  attempted  for  the  converging- 
diverging  nozzle  example.  In  one,  two  interior  stations  were  isolated,  and 
the  forward  divided  difference  was  used  on  all  variables  that  were  differentiated 
with  respect  to  x.  Letting  two  consecutive  p  values  be  called  R1  and  R2, 
two  consecutive  m  values  be  XM1  and  XM2,  and  two  consecutive  e  values 
be  El  and  E2,  a  system  of  six  state  variables  would  result.  However, 
temperature,  which  could  be  considered  constant  but  is  actually  a  slowly 
varying  function  of  e,  m,  and  p,  and  d(ln  A)/dx,  a  constant,  must  be 
accounted  for  at  the  two  points.  By  setting  their  derivatives  with  respect 
to  t  to  zero,  these  constants  can  be  input  to  STAN  along  with  the  state 


variables.  Let  XK1  and  XK2  be  consecutive  values  of  d(ln  A)/dx  and  TMP1 


and  TMP2  be  temperature  values,  and  denote  the  derivative  of  a  variable 
Z  by  Z*,  then  the  resulting  input  to  SERIES  is: 


R1 . =- (XM2-XM1) /DX-XM1*XK1 ; 

R2 . =- (XM2-XM1) /DX-XM2*XK2 ; 

XM1. =-XMl*XMl/Rl*XKl+( (XM1/R1) **2-R*TMPl) * (R2-R1) /DX 
-2.*XM1/R1*(XM2-XM1)/DX; 

XM2 . =-XM2*XM2/R2*XK2+( (XM2/R2) **2-R*TMP2) *(R2-R1) /DX 
-2.*(XM2*XM2/R2)*(XM2-XM1) /DX; 

El . =-XMl* (El /R1+R*TMP1) *XK1+(XM1*E1) /Rl**2* (R2-R1) /DX 
-(E1/R1+R*TMP1)*(XM2-XM1)/DX-XM1/R1*(E2-E1)/DX; 

E2 . =-XM2* (E2/R2+R*TMP2) *XK2+XM2*E2 /R2**2* (R2-R1) /DX 
- (E2/R2+R*TMP2) *(XM2-XM1) /DX-XM2/R2*(E2-E1) /DX; 

TMP1.=0. ; 

TMP2.=0. ; 

XK1.=0. ; 

XK2 . =0 . ; ; 

where  R  and  DX  are  constants  and  the  known  equilibrium  values  from  a 

_  2 

test  run  can  be  read  in.  Since  TMP  =  (e  -  m  / 2p)/pc ,  TMP1.  and  TMP2.  can 
also  be  entered  by  differentiating  this  expression,  but  results  will  be 
similar.  Appendix  \  contains  the  output  of  SERIES  for  this  input. 

In  order  to  avoid  using  the  same  derivative  with  respect  to  x, 
a  system  based  on  one  x  point  with  constant  input  partial  derivatives  was 
also  tried.  This  system  uses  the  variables:  RX  for  p;  XM  for  m;  E  for  e; 
RDX  for  3p/3x;  XMDX  for  3m/3x;  and  EDX  for  3e/3x;  TMP  for  temperature,  and 
XK  for  d(ln  A)/dx. 
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RX.=-XM*XK-XMDX; 

XM. =-XM*XM/RX*XK- (R*TMP- (XM/RX) **2) *RXDX 

-2 . *XM/RX*XMDX; 

E . =-XM* ( E / RX+R*TMP ) *XK+XM*  E / RX*  *  2  *RXDX 

- (E/RX+R*TMP) *XMDX-XM/R*EDX ; 

RXDX.=0. ; 

XMDX.=0. ; 

EDX.=0. ; 

TMP.=0. ; 

XK.=0. ;; 

Appendix  B  is  the  output  of  SERIES  for  this  input. 

Both  systems  were  tested  against  the  explicit  Euler  solution  of  the 
corresponding  initial  value  problem  (SERIES  also  generates  a  FORTRAN 
subroutine  DIFFUN(T,  Y,  DY)  which  fills  the  array  DY  with  the  derivative 
evaluated  at  u(t)  where  t  =  T,  u  is  in  array  Y) .  It  was  discovered  that 
the  radius  of  convergence  of  the  power  series  contracted  sharply  for 
values  past  the  throat  of  the  nozzle,  so  only  values  between  the  inlet 
and  the  throat  can  be  analyzed  using  STAN.  Table  I  gives  initial  values 
that  were  picked  for  analysis.  Note  that  the  throat  is  at  x  =  .26 
where  d(ln  A)/dx  =  0. 

Both  of  these  systems  were  run  with  both  sets  of  initial  data,  and 

the  resulting  contractivity  regions  are  displayed  in  Figure  4 .  These  were 

-12 

only  achieved  for  At  of  5*10  ,  and  do  not  correspond  to  the  expectations 

of  results  from  test  runs.  Also,  they  are  identical  for  both  the  analytic 
and  numerical  solution,  which  suggests  they  are  actually  an  artifact  of 
the  program  STAN.  This  can  be  seen  to  be  the  case  since  the  first  step 
of  generating  the  stability  region  about  an  equilibrium  point  U*  is  to 
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Figure  4.  # Non-linear  stability  analysis  of  the  nozzle 
problem.  Solid  line:  m  vs.  e  with  center  at 
(16.5,  355587.).  Dashed  line:  p  vs.  m  with 
center  at  (0.,  16.5).  Both  axes  of  length  1 
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3.  Research  Findings 

Relative  to  the  specific  research  outlined  in  the  original 
proposal  on  this  research,  the  following  findings  are  of  interest. 

1.  A  study  of  recent  literature  on  finite  elements  [18], [19] 
reveals  that  use  of  finite  element  techniques  are  not  easily  adaptable 

to  nonlinear  systems  in  several  dimensions,  especially  when  the  system  is 
designed  to  be  easily  changed  as  different  reentry  vehicle  configurations 
are  tested.  The  analysis  of  minimum  and  maximum  step  size  found  in  [12] 
is  really  inadequate  in  the  finite  difference  case,  depending  as  it 
does  on  linearization  of  the  differencial  system,  and  a  similar  analysis 
in  the  finite  element  case  seemed  both  beyond  the  scope  of  the  intended 
research  and  not  very  fruitful. 

2.  Since  the  algebraic  amplification  matrices  involved  in  the  linear 
stability  analysis  of  the  methods  under  investigation  did  not  point  out 
the  experimentally  observed  instability,  it  was  inappropriate  to  develop 

a  complicated  algebraic  manipulation  package  to  compute  such  matrices. 

3.  After  working  with  the  two  model  packages  HVSL  and  the  converging- 
diverging  nozzle,  it  was  concluded  that  the  form  in  which  the  updates  to 
the  dependent  variable  vector  is  derived  (the  "delta"  form) ,  resulting 

in  parts  of  the  numerical  method  being  computed  at  various  stages  and  places 
in  the  program,  would  make  it  difficult  to  submit  these  codes  to  use  by 
standard  packages,  [20],  [21]  which  usually  require  a  single  subroutine  such 
as  DIFFUN(T,  Y,  DY)  to  compute  the  first  derivative  array  DY  given  the 
state  variable  array  Y  and  the  independent  variable,  T.  On  the  other 
hand,  the  "delta  form"  can  be  easily  adapted  to  most  implicit  multi-step 
correctors 

k  k 

U(x  )  =  E  o  U(x.  .)  +  Ax  E  B  U(x,  .), 

J  i=l  1  -1-1  i=0  1  -'_1 
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search  from  U*  along  a  particular  line  in  (u^,  u^)  space  emanating 

from  U*  searching  for  the  initial  conditions  of  the  bisection  method,  i.e., 

find  i  such  that  AV(U*  +  2*(e^,  e2))  has  different  sign  from 

AV(U*  +  2i+^(e^,  e^)),  where  (e^,  62)  represents  the  unit  vector  in  the 

search  direction  for  the  two  variables  being  changed,  u^  and  u2*  The 

variable  i  varies  from  -3  to  4.  In  most  applications  of  ordinary 

differential  equation  stability  regions,  the  various  state  equations  are 

not  highly  coupled,  and  when  two  equations  are  coupled,  one  usually 

attempts  to  find  the  stability  region  using  these  two  variables.  For 

discretized  partial  differential  equations,  however,  the  variables  are 

necessarily  highly  coupled,  and  apparently  a  change  on  the  order  of 

100  percent  creates  an  immediate  overflow.  Thus,  the  actual  stability 

-12 

region  for  At  =  5*10  occurs  because  AV  =  0  when  At  is  made  so  small 
that  VCu^O)  ,u2(0))  =  V(u^(At)  ,u2(At))  by  either  numerical  or  analytical 
techniques. 

To  test  this  hypothesis,  one  value,  at  x  =  .24,  of  p,  of  m,  and  of 
were  each  changed  to  10  times  their  original  value  and  the  original, 
linearized,  converging-diverging  nozzle  program  was  run.  In  each  of  the 
three  cases,  overflow  stopped  the  calculation  by  the  third  time  step. 
Therefore,  STAN  must  be  modified  to  take  into  account  the  coupling  of  the 
state  variables  when  analyzing  the  stability  of  systems  of  partial 
differential  equations.  This  modification  is  currently  being  made.  If 
successful,  it  would  open  the  way  to  automatic  stability  analysis  of 
complicated  systems  of  nonlinear  equations  to  allow  researchers  to  choose 
which  numerical  method  is  most  appropriate,  from  a  stability  standpoint,  to 
integrate  their  parabolic-hyperbolic  discretized  system. 
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with  only  the  additional  storage  for  carrying  the  U(Xj  ^) ,  U^(Xj 
terms.  Also,  methods  for  changing  step  size  Ax  and  even  changing  from 
one  formula  to  another  could  be  adapted  from  the  packages.  However, 
since  the  finite  differences  are  usually  only  first  order  accurate  the 
effort  on  nonlinear  systems  of  a  high  order  method  in  the  independent 
variable  versus  a  low  order  method  in  the  spatial  discretization  is 
not  yet  understood.  Certainly  in  the  linear  case  there  is  a  stability 
argument  against  it. 

4.  The  best  tool  now  available  for  the  analysis  of  a  nonlinear 
parabolic  system  could  be  STAN,  provided  an  understanding  of  how  to 
reduce  a  discretized  system  to  one  of  a  convenient  size  for  such  an 
analysis.  Certainly  the  cost  of  both  computer  time  and  programmer  time 
of  entering  the  entire  discretized  system  is  prohibitive,  yet  the  two 
attempts  to  enter  significant  subsystems  did  not  yield  sufficient  information 
to  show  the  value  of  this  analysis  technique. 

However,  experimental  evidence  exists  that  the  numerical  methods 
currently  used  to  model  high  speed  flow  on  a  cone  is,  at  best,  marginally 
stable,  and  the  results  suspect.  Continued  research  should  be  undertaken 
to  provide  an  understandable  method  for  directly  analyzing  the  stability 
properties  of  parabolic-hyperbolic  systems,  and  comparing  them  to  the 
stability  behavior  of  numerically  generated  solutions,  and  to  chose,  when 
appropriate,  more  accurate  numerical  methods  that  do  not  require  significantly 
larger  storage. 
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appendix  a 


Output  of  SERIES  for  divided  difference  formulation 


SUBROUTINE  01 FFUN ( T  *Y , DY ) 

0 1  MENS  I  ON  QY ( 20 ) »  Y ( 20 1 

_ DA TA  R/1 7I6./.XKI /7.867 1 1  /_,  XK  2 /_7_*__5? 839/tC V/ 4  2  9  0_.  / ,_DX /  . 01/ 

DY( 1) =-( YI4)-Y(3) )/DX-Y< 3)*XKl 
DY( 2) =-(Y (4>-Y(3) ) /DX-Y ( 4 ) ' XK 2 
iJY  (  3  )  =-Y  I  3  )  *  Y  (  3  )  /  Y  (  I  Y  /YJ_L  >  > **2-R  » Y  (  7  )  ) »  (  Y  (  2 )  -Y  (  1 )  )  /DX- 

♦  2 • *  Y ( 3 ) / Y ( 1 ) * (Y 1 4 ) — Y ( 3 ) ) /DX 

DY( 4) =-Y(4) *Y(4)/Y(2)*XK2*( ( Y(4»/Y( 2) ) **2-R*Y ( 8) ) * ( Y  ( 2 1  — Y (111 /DX- 

+  2.MY(4)/Y(2 ) )* ( Y(4)-Y(3i ) / 0  X_  _ _ _ 

0Y<  5)=~Y(3I*(Y(5I  /Y(  U+R*Y<  7)  |*XK  U  {  Y  (  3  f*Y  (  5  H  /  Y(I)  **  2*  (  Y  (  2  l-Y(  1 1 

♦  ) /OX-1 Y(5)/Y( 1)+R*Y(7| )*  <  Y ( 4  I  — Y  C  3  > ) /DX— Y I  3  > /Y  C  1  )* ( Y ( 6 ) -Y ( 5 ) ) /OX 
0Y(  6)=-Y(  4)«(  Y  (6  )_/Y  (  2  )_+R  'Y  (  8 )  )  *  XK 2* Y  l  4  )_*  Y  (  6) /Y(  2)**2 *  I  Y  I  2  >  —  Y  <  1 1  >  / 

♦  OX*  (  Y  (  6  )  /  Y(  2  >  +  R*Y  (  8  ) )  *  (  Y(4)"-Y  ( 3  )  )  /OX-Y(  4)~/Y  (2)*(Y(6)-Y(  5)  I  /OX  ' 

0Y<  7)  =!  U-Y(3)*(Y{  5)/Y(  1 )  ♦R*Y  ( 7  ))  *XKI  +  (  Y(  31*  Y(  5)  »/Y(l)**2*4  Y 1 2  » — Y 

♦  11)  )/DX-(  Y(  5)/Y(l  )+K«Y17  J)*m  4)-Y(  3)  l/OX-Y  C  3  > /Y  (  1)«(Y(6)-Y(  5)  l/OX 

+  )-Y(3)M-Y(  3)  *Y(3 )/Y(  l)*XKl  +  ( (Y( 3) /Y( l)  »**2-R*Y( 7) )  *(  Y(  2 1 -YU)  ) /DX 
+—2 .  -rY  (  3  )  /  Y(  1)  M  Y14J-YI  3»)  /OXI/Yl  1)  )/YC  l)-<  Y(  5 ) -Y  (  3 )  *•  2/ Y  (  1 )  )*<-(  Y( 
♦4)-Y( 311 /DX— Y (3)*XKll/Y(l)**2)/CV  _ 

DY  (  8  j'=  (  (  (  -Y  (  4  )  «  (  Y  (  6  1  /Vi  2  )  f  R*Y  {  8  )  )*XK2+Y  ( 4) * Y( 6 ) /Y ( 2  >**2 * ( Y( 2 )-Y( 1 
+  ) ) /0X-( Y( 6)/Y<  2)+R«Y{ 8 >) * ( Y < 4 » -Y ( 3)) /DX-Y ( 4 ) / Y ( 2 ) * ( Y ( 6 ) -Y { 5 )) / DX ) - 
+  YC4)  M-Y(4)*Y(4  )/Y(  2)*XK2+(  ( Y  (41 /  Y(  2) )  **2-R«Y(  8D*  (  Y(  2)-Y<  11 )  /OX-2 
+  .M  Y(  4)/Y(  2)  )M  Y<  4)  — Y  (  3)1  /DX) /Y  (2  )  )  /Y < 2 ) -( Y ( 6 1 -Y ( 4 ) ** 2/ Y(2 ) ) * ( - ( Y < 
♦4 ) — Y 1 3) )/QX— YI4)*XK2)/Y12)**2 )/CV 

OY ( 9 ) =  1 . _ 

RETURN 

END 


SUBROUTINE  SOL C T. YO, YNEW » INOI 

DIMENSION  YO ( 20 ) »  YNEW(  20),ZZZB(  20  ) » DTP  AKE  (  20)  ,  DR  I  (  201  ,Rl<20), XM2 

_ ±f  20),XMl(20).TST(20)t  TSU(20  )  ,_TSVJ  20 J  , TS W( 20 ) . TSX { 20 )  ,DR2(  20). R2 

+  ( 20) , TS 1 ( 20) ,TS2( 20) . QXMl (201 ,TS3< 20) ,TS4C 20) , TS5( 20) , TS6I 20) , TS7 
+(20),TS8(20I,TMP1(20),TS9(20) ,TTS( 20) ,TTT( 20) ,TTU( 20) , TTV( 20) »TTW 

_ +<_2  0),TT  XJ  2  0)  ,  TTY(  20) .  TTO  (20  )  »  TT  l  (  20 )_,  TX?i2  0 )  .  DXM2<  20)  ,_T  T _3J  2  0  L»  TT4 

+(20)»TT5(20)»TT6(  20) , TT 7( 20 ) , TT 8( 20 ) , TMP2 ( 20 J , TT9( 20) , TUS ( 20) » TUU 
+  (  20 )  »  TUV(  20) »  TUW<  2  0)  *  TUY(  20) »  TUO(  201  *  TUK  20)  .  TU2  (  20)  »  OE I  (  201  »  El 
+  (  20)  ,TU3(20)  ±  TUf>JL2  0 ) , TU6 ( 2  0) , TU 7 ( 20 ) ,  TUB  (JO )  ,  T  U  9(20  ).TVS(  20)  >TVT 
♦(20) ,TVV(20),TVW( 20) ,TVX( 20) , TV2( 20) , TV3(  20) , TV4< 20) * E2 (20 > ♦ TV6 
+<20>,TV7(20),TV8( 20) , T V9( 20 ) , 0E2( 20) . TWS( 20) . TWUI 20) , TWV(  20) , TWW 

_ +( 20 ) ,  TWX(  20 ) »  TW Y (  2 0 )  ,  TWZ  (  2_0_»  »  T WO^t  20 )  ,  TW2X20 )  ,TW3(  20)  .  TW4(  20)  ,  TW9 

+  (20) ,TXS( 20).fXT( 2  0) »TXW (20). TXX ( 20 ) » TXY ( 20) , TZ6(20 ) , TZ 7( 20) ,TZ8 
+  (20) , TZ  9( 20) . TOS( 20) , TOT ( 2 0) , TOU ( 20) , T00( 20) , T02 ( 20) , T03( 20 ) , T04 
+  (  20) ,T3U( 20 ) . T3V(  20) . T 3W( 20 ) . T3X( 20 ) ♦ T3Y ( 20) , T3Z ( 20) , T30( 20) .T36 
+(20),T38(20),T39( 20 ) , T 4 S ( 20 ) . DTMP1 ( 20 ) , QTMP2 ( 20 ) , T FAKEl 20 ) 

DATA  R/ 1716./ i  XK1  / 7. 3 6 7 1 1/ , XK 2/ 7 . 59839/ , C V/42 90 . / , DX/ . 0 1/ 

E  PS=  1 . 0E-6 _ _ _ _ _ 

R  f(  1 )  =  YU (  I) 

R2( 1) =Y0(2) 

XM1  (  I )  =YU  (  3 ) _ _ 

XM2 ( 1 ) =Y0( 4 ) 

E 1 ( 1) =Y0( 5) 

_ E2  (  I )  =  YO  (  6 ) _ _ 

TMP  1  ( 1 ) =Y0( 7) 

T  M  P  2 ( 1 ) =YQ( 8) 

T  FAK  E ( 1 )=Y0<9) 

I  ND  =  0 

DO  1  111=1,19 
N 1 1 1  =  I  I  I 
IIIl=lIl  -  1 

T  ST  C I  I  I )  =  XM2  (  1 1 IJ-XMK  III) 

TSUI  III  )=TST(III)/DX 

TSVTl  I  I  )  =  -TSUn  II  )  ' 

T  SM ( I  II )=XMl( II I)+XK1 
T  SX  (  I  I  I )=TSV( II I)-TSW( III) 

DRl(IirT=T$X(JII) 

RHIII  ♦  l)  =DR1(  1  II  )/FLOAT(  II  I) 

TSl ( I  1 1 ) =XM2( I  I  I) «XK2 
TS2(  1 1 1 )  =  T  SV(  Iin-TSHI  II ) 

DR2 ( I  I  I )  =  TS2( III) 

R  2  < III  +  1) =DR2( I  1 1 )/FLOAT(  III ) 

T  S3 ( I  I  I ) =0. 

DO  100  JJ J=l, I  I  I 

100  TS3  (  I  1 1  )  =TS3(  1 1  I)  +XM1(  JJJ  )  +  XMK  I  I I-JJJ  +  1) 

IF  (III. EQ. 1)  GO  TO  101 

TS4( I  II )=TS3(  1 1 1 ) — TS4 ( l)*Ri( I  II I 
IF ( II I.EQ.21G0  TO  102 
00  103  J J J=2» I  I  II 

103  TS4( I  I  I >  =  TS4( I  I I)-TS4( JJJ)*R1( I  I I-JJJ  +  1) 

102  TS4( I II)  =  TS4( 1 1 1 ) /Rl  ( l ) 

G  O  TG  '  104'  . . 

101  TS4(  I  II )  =  TS3( 1 1  I) /Rl(  1 1  I) 

104  CONTINUE 

TS!>(  I  I  I  )  =  TS4(  H  I)  «XKl 
TS6(III)=-TS5(III) 

IF  (III.EQ.l)  GO  TO  105  _ 

TS7< 1 1 1 )  =  XM1 ( 1 1 I)-TS7( l)»Rl( I  II  ) 


IK  I  I  I .EQ.2  )G0  TO  106 
00  107  J J Js2* 1 1 1 1 

107  _ TS7J  M  IJ^JS7_jn_n-TS7(  JJJ|*R1I1 1  L“JJ«J  +  IJ_ 

106  T  S7<  I  I  I)=TS7I  1 1  I)  /RKl) 

GO  TO  108 

105  _ TS7 1 II I  )  =XM1 1  IIII/RK  IIII _  __  _ 

108  CONTINUE 
TS8(III)=0. 

_ 00  109  JJJ=1, 1 1  I _  _  _  _ 

109  T  S3 1  III  >  =  TS8(  III)+TS7(  JJ J > *TS7( 1 1 l-JJ J+ 1) 
TS9( 1 1  I )  =  T  M  P 1 ( I  II  )*R 

_ TTS( II I )  -  T  S  8  (  IID-TS9 (III) _ 

TTTt I  II )  =  R2(I I  I )— R 1 ( I  I  I) 

TTUI 1 1 1 )  —  0 • 

_ oo  no  jjj=i tin _ 

110  TTUI  II  I)  =  TTUI  II  D+TTSI  JJJI*TTTl  I II-JJJ+1) 
TTV(  III)  =TTU (  IID/DX 

TTW( I 1 1 )  =  TS6l III)+TTVt III) 

TTX ( 1 1 1 )  =  XM1( 1 1 1 ) *2 . 

IF  ( I1I.6C.I)  GO  TO  111 

TTY l I  1 1 )  =  TTX( I  I  I) -TTY( 1)*R1(I II) 

IF  C II I.EQ.2IG0  TO  112 
DO  113  JJJ=2,III1 

113  TTY  I  II [)  =  TTY( II  I ) -TTY ( J J J ) * P l (II I-J J J* 1 ) 

112"  TTY  (  I  I  I )  =  TTY  (  1 1  I )  /Rl  f  l ) 

GO  TO  114 

111  TTY(  I  I  I )  =  TTX( 1 1  I) /Rl(  I  1 1) 

114  CONTINUE 
TTO( I  1 1 ) =0. 

DO  115  JJJ=1,III 

115  TT0( III)=TTO( III)+TTY( JJJ)>TST( III-JJJ+1) 
TT 1 ( III J-TTOC III) /OX 

TT2( 1 1 1 )  =  TTW( 1 1 1)-TTl( I  I  I  ) 

DTX  MK  1 1  17  =  T  T  2  ( 1 1  rr~ 

XMKIII  +  l)=OXMi  (  III)/FLOAT(  III) 

TT3 ( 1 1 1 ) =0. 

DO  116  JJJ=l, m 

116  TT3 ( I  I  I  I  =  TT3  < III) +XM2( JJJ)«XM2( I  I I-JJJ+1) 
IF  ( III.EQ.l)  GO  TO  117 

T T4 (TTH  =  TT 3 TTI  lT*TT4"("l  J *R2 rfTTT 
IF ( 1 1  I .EQ.2 )G0  TO  118 
DO  119  J J J=2t  II  II 

1 19  TT4(  1 1  I  )  =  TT4  ( I II)  -TT4Tj  JJ  )  «  RTHTI-J  J  J  +  l  ) 
118  TT4( IIl)=TT4( 1 1 1 ) /R2 ( l) 

GC  TO  120 

117  T  T 4TTni  =  T T  3TTI  I )  / R2  (III  1  “ 

120  CONTINUE 

T  T  5 ( I  1 1 ) =TT4 ( III) *XK2 

t  t  6TTitt=  -  t  t  5  rnr  > "  ~ 

IF  ( III.EQ.l)  GC  TO  121 
T T 7 (  I  I  I )  =  XM2(  I I  I )  "TT7(  l )  "*R2  (III) 

I F ( 1 1  I . E  Q  .2 I  GU  TO  122 
00  123  JJJ*2fIIU 

123  T  T  7  (  I  IJ  » =  TT_7J  1 1 1)-TT7(  JJJ)*R2(I  I  I-JJJ  +  l  I 

122  f  t 7 (  1 1 1 ) =  T T 7 (  1 1 1 ) /R2  ( T) 

GC  TO  124 

121  _ T  T  7  (  II  I )=XM2( IT  I) /R2( 1 1  I) 

124  CONTINUE 


1 


TT8(III)>0. 

00  125  JJJ=l,IIl 

12 5  _ T T a^I  I  I  )  =JX8 ( 1 1 1  > »TT 7(JJJ>*  TT7(  I  I l-J JJ+l ) 

T T9 ( III) *TMP2 (ill )*R 
TUS( I II  I  =  TT 8 < III)— TT9I 111) 

_ TUUI  I  1  I  )  =  0 . _ _ _  _ _ 

00  126  JJJ=l, II  I 

126  TUUUlU=TUUt  IIIH-TUSI  JJJ)*TTT(  IIl-JJJ+1) 

tuvi  1 1  ii  =  xu  uiiii)/  ox _ _ 

Tuv»(irn=Tf6(iii)+fuv(iin 

TIJYUII  )  =  TT7<  1111*2. 

TUO ( I  I  I ) =0. 


00  127  JJJ=1, I  I  I 

12  7  TUOI I  II »  =  TUOIII I )  +  TUY(  JJJ)*TST(  I  I I-JJJ  +  1 ) 

_ TUI ( I  1 1 ) =TU01  1 1  I) /OX _ 

TU2( 1 1 1 )  =  TUto ( 1 1  I) -TUI  I  1 1 1 ) 

0  X  M  2 (  III) =TU2(  III) 

XM2  (  I  1 1  +  1)=DXH2(  III) /FLQATC  III)  _  _ 

'  IF"  (  I  II  .EQ.l)  GO  TO"  128  ”  . " 

TU3<III)  =  E1(III)-TU3(1I*R1(  III) 

_ IF<  111. EQ. 2) GO  TO  129 _ 

00  130  JJ J=2, 1 1  II 

130  TU3  (  I  1 1 )  =  TU3 (  I1D-TU31  JJ J ) *Rl ( 1 1 1-J JJ  +  l ) 

129  TU3  (  I  II )  =TU3(  I  ID  /Rl(  1  ) _ 

GO  TO  131 

128  TU3( 1 1  I ) -El  II  1 1  ) / R1 ( I  I  I ) 

131  CONTINUE 

TU5 ( I  II)  =  TU3( 1 1 1 ) ♦tS9( II 1 1 
TU6( I  I  I ) =0 . 

00  132  JJ J»l, 1 1  I 

132  T  U6  <  I  1 1  )  =  TU6  <  1 1  U  ♦  XMl  <  J  j  J  J  *  f U  5  <  1 1 1  -  J  J  j  + 1 } 
TU7 1  I  1 1 ) =TU6I II  I) *XK1 

TU8 ( 1 1 1 )=-TU7( III) 

TU9 ( 1 1 1 ) =  0 • 

DO  133  JJJ=1,  III 

133  TU9I  III  )  =  T'J9<  II  I)  ♦XMK  JJJ)*E1(  I  I  I-JJJ  +  1) 

T VS ( I 1 1 J  =0 .  . . 

00  134  JJJ=1,III 

134  TVS!  1  1 1 )  =  TVS  C  III)  +  R1I  JJJMRll  II  I-JJJ+1) 

IF  (III. EQ.l)  GC  TO  135 

TVTI I 1 1 )=TU9( 1 1 II-TVTI 1 )*  TVS! 1 1  I ) 

I F ( II I.EQ.2IG0  TO  136 
DC  137  JJJ=2, III1 

137  TV T ( I  1 1 )  =  TVT l I  I  I ) -TVT ( J J J ) *TVS( I  I I-JJJ+1) 
136  TVT( III l=TVT( I  I  I) /TVS( 1) 

GO  TC  138 

135  TVTI  1 1 1 )  =  TU9( II I) / TVSIII I  ) 

138  CONTINUE 

T VV ( I  I  I ) =0. 

DO  139  JJJ  =  1,  1 1  I 

139  T VV ( I II)  =  TVVXIII) +  TVT<  JJJ)“TTT( III-JJJ+II 
T  Vim  (  1 1 1  )  *  TV  V I  1 1 1 1  /O  X 

T VX ( I  I  I ) =TU8I  I  I  I) +  TVW<  III) 

_ TV2( I 1 1 )=0. _ 

JO  140  JJJ=l,  I  I  I 

140  TV2( I  1 1 )  =  TV2l 1 1  I) +  TU5  ( J J J ) *  T  ST (  I  I I-JJJ+1) 
T V 3 1  III) =TV2( II I) /OX 

1V4'(  III  )  =  TV  X  (  I  I  II-TV3I  III) 


141 


Ji... 


TV6Uin=E2nii)-Ei(iin 
TV7(  nn=o. 

OQ  141  jjj=i,ih _ _  _  _ 

T  V  7  (  I  1 1 )  =TV7< I  II) +TS7I JJJ)*TV6I II I-JJJ+1) 
TV8( I  I I)=TV7( I  1 1) /OX 

-LV'LL! li  I^TV4(_I  IJ)-TV8(l  II I _ _ _ 

0  E 1 ( II l)=TV9(  III) 

e i (  iii  +  i)=dehiii)/float(iii) 

IF  I  III. EG. 1)  GC  TO  142 _ _ _ 

TWSI  I  II)=62(  1 1 1  )-TWSm*R2I  III) 

I F t  1 1  I.EQ.2)CiO  TO  143 
_ 00  144  JJ  J  =  2,  I  I  1 1 _ _  _ 

144  TWSU  m=TnS<  1 1 1)  -TwSI  JJJ)«R2<  fl  I-JJJ  +  1  > 
143  TwS(III)=TWS( III)/R2(1) 

_ GC  TO  145 _ 

142  TwS{IlI)=E2(Iin/R2(III) 

145  CONTINUE 

_ TWU(  1II)  =  TWS(  II  I) +TT9(  Ill  ) _ _  _ 

r^v'ii  1 1 )  =  o  Z 

DO  146  JJ J=1 ,  I  I  I 

146  TW  V  (  I  1 1  )  =  TW V I  I  I  1) +  XM2 ( J J J )  “ TUU(  I  I 1-JJJ+l) 
TWWII 1 1 )=TWV<  1 1  2) *XK2 

t  smx  1 1 1 1  i=-mwi  iii  ) 

_ T WY ( I  I  I ) =0. _ 

DO  147  JJJ=1,  I  I  l 

147  TWY  (  1 1 1  )  =  TWY {  1 1 1) +XM2 ( JJJ)*E2II I l-JJJ  +  1) 

_ TUZII I  I ) =0  . _ 

00  148  JJJ=1, III 

140  TWZ(III)=TWZ<  IID+R2I  JJJ)*R2(  II  I-JJJ+1) 

IF  ( III. EG. 1)  GC  TO  149 

7wo<iii)  =  7wy(  iin-f*on)*~Trtinu) 

I F (  1 1 1 • £0.2 )GU  TO  150 
00  151  JJJ=2»III1 

151  TwOll  1 1  )  =  T  WOT  1 1  I)  -TWOUJJ)  +TWZTI  II-JJJ  +  1) 
150  TwO( I  II )=TWO< I  II) /TWZI 1) 

GO  TO  152 

1 19  TWO< I II )»tWY( III) /tWZI III) 

152  CONTINUE 
TW2 I  1113=0. 

DO  153  JJJ=1,III 

153  TW2I  III  )  =  TW2( II I)+TWO( JJJ)*TTT( 1 1 1-JJJ+l) 
TW3 I  I  1 1 )=TW2(  1 1  1) /DX 
TW4 ( I  1 1 )  =  TWX ( I  II) +TW3II  IT) 

TW9I  I l I )  =  0 . 

00  154  JJJ=l,III 

TW9(  III)  =  TW9  f  III )  +  TWU  (  J  J  J  )  *  T  STTlI  I-JJJ+1  f 
TXS ( I  1 I )  =  TW9( 1 1  I) /OX 
TXTI I  1 1 ) =TW4( I  I  I ) -TXS ( III) 

tx¥(Th  )=o. . 

DO  155  JJJ=l, I  I  I 

T  X  W  (_I  1 1 1  =TXW(  1 1 1 )  +TT7(  J  J  J  )  *  T  V6  {  1 II-JJJ+I) 

T XX  (  1  H  )  =  TXW(T  H)  /OX 
TXY(III)«rXT( II I)-TXX( III) 

0E2  (  in  )  =  TXY(  1 1  I) 

E  2  (  III  ♦  l)=0E2( III )  /  FLOAT! II  I) 

T Z 6 (  I  I  l  ) =0 . 

00  156  JJ J=lf I  I  I 

T  Z  6  <  n  I  j»TZ6(  I  1 I) +XM1I  JJ  J)  mTT2(  I  l  I-JJJ+ll 


154 


_ 15  5 

I 


1  66 


IF  { III .EG. 1)  GC  TO  157 

TZ7(  I  II  }=TZ6t  III)-TZ7m*Rl(I  II) 

_ I F  <  1 1  1  .EQ.2  )GO  TO  158 _ _ 

bG  159  J  j  J  =  2  *  I  I  II 

159  TZ  7  (  I  I  H  =  TZ7<  I  1 1 ) -TZ 7 < JJJ 1 « R1  (II  I-JJJ+l ) 

158 11 7 ( I  I  I ) =  TZ 7 1  I  I  I) /R 1 U  ) _ 

GC  TO  150 

157  T  Z  7 ( I  I I)  =  TZ6( I  I  I) /Rl< 1 1  II 

ij»q _ coNTj_Nye_ _ _ __ _ 

fz8(lni  =  TV9(  iiii-f'z7uin 

IF  (III .EJ.l J  GC  TO  161 

TZ9mi )  =  TZ 8(  I  1 1 1 -TZ9( 1 )  "R1 ( I  II  ) 

IP<II I.EQ.2 1G0  TO  162 
DO  163  J J  J  =  2  f  1 1  1 1 

1  <>3  TZ9( II  I ) =TZ9t I  I  I) -TZ9( JJJ)*R1  (  I  I  I-JJJ+l 1 
162  TZ9( I  1 1 )  =  TZ9(1 1 1) /Rill) 

GC  TC  164 

1 d 1  T Z  9(1 II )  =  TZ 8(1111/^1(1 m 

164  CONTINUE 

T OS ( I  I  I  1=0. 

DC  165  JJJ=1,  II  I 
155  T  OS  (  I  1 1  I  =  TO  S  (  1 1 11  +XMU  JJJI^XMII  1 1 1- j  JJ  +  1 1 

IF  (III.EQ.il  GC  TO  166 
TOT ( I  1 1 )  =  TO S (  1 1 I)-TOT( 1)*R1(  I  II ) 

I  F(  II  I.  EQ.2  I  GO  TO  167' 

OC  168  JJ J=2,  1 1  II 

168  T  OT ( II I  I =  TO  T (  1 1  I ) -TOT  I J J J ) *R1 (I  I I-J JJ  +  1 1 

167  roT(iin  =  TOT<iin/RnV> 

GC  TO  169 

166  TOT( I  I  I )  =  T0S<  I  I  l) /Rl( I  III 


T 02  (1111  =  TOO (  I  IT)  /T  VS  Cl  I  I I 
CONTINUE 

T03( I  I  I  I =TZ9( III I-T021 1 1 1  ) 

T04 ( I  1 1  I =T03 1  I  I  I) /CV . . 

DTMPK  III  I  =  T04 (  II  I  1 

TMPKIII  +  11=DTMP1(  IIIl/FLOATI  III) 

T3U(  I  1 1  1  =  0. 

DO  175  JJ J=l,  II  I 

T3U(I 1 1 )  =  T3U( II II +XM2 ( J J J 1  *  TU2(  1 1 1-JJJ+l) 

IF  (  HI.Ei3.il  GO  TO  176 

T  3  V ( l  1 1  I  =  T3U ( I  I  II -T3V(  l)*R2( I  II  I 

IF(  1 1  I. E  J. 2  1  GO  TO  177 

DO  1 7  8  J J j=2* I  I  II 

T  JVIII I  )  =  T3V(  I  I  U-T3V(  JJJ  >*R2(  l  I  l-J  JJ  +  1) 

T 3 V ( 1 II) =T3V(  II  I) /R2(  1) 

GO  Id  179 


1  71 
174 


175 


l  78 
1  77 


19  2  T3  8C I  1 1 )=T38( 1 1 1 ) — T 38 ( JJJ)*TWZ( 1 II-JJJ+1) 
19l  T  3  8  I  1  I  I  )  =  T  3  8  f  1 1  IJ /TwTT 1  ) 

GO  TO  193 

L90  T38(  I  I  I  )  *  T3  6  (  I  I  I)  /  TinlZ(  III) 

193  CONTINUE 

T 39(11 I)  =  T3X( II I)-T38( 1 1 1 » 


195 


196 


ZZZZ  1 
I F  (  J  J 
Jill  2 
CGNTI 
ZZZZ1 
I  F  (  Z  Z 


=  ZZZ 

J.LT 

^ZZZ 

MUE 

=  E  PS 

ZZ2. 


Z1  +R2IJJJ) 

.1  1 1 -4 )G0  TO  195 
Z2  +  ABS ( R2 ( J J J  ) ) 


*  (  ABS l ZZ  ZZ  1 )  + 
GT. ZZZZI ) GO  TO 


) 


ZZZZ1 
Lilli 
on  19 


=o. 

=o. 

6  JJ 


zzzzi 

IF(  JJ 

lllj.2 
CGNTI 
ZZZZI 
IF(  ZZ 


=in 

J.LT 

=  lll 


J=1^I I  II  _ _ 

Z1  ♦  XMilj'jJ) 

.1  II-4  )G0  TO  196 
12  +  ABS ( X»1 ( J JJ  )  ) 


NOE 
=  EPS 
111. 


MABSIZZZZll  + 
GT. ZZZZI )G0  TO 


> 


till  1  =  0. 

ZZZZ2=0. 

_ 00  197  JJ  J=lt  I  I  II 

11111=11111  +XM21JJJ) 
IFUJJ.LT. I  1I-4JG0  TO  197 

_ ZZ  1 12 -  III  11  »  A3S(XM2(JJJ)> 

197  CGNTI HUE 

ZZZZ1  =  EPS;*(ABS(ZZZZ1)  ♦  l.) 

_ IF(  ZZZZ2.GT.ZZZZHGQ  TO  1 

ZZZZ1=0. 

Ill  12-0 . 

DO  198  JJJ=lrIIIl 
ZZZZ1  =  ZZZZ1  +  6KJJJ) 

I F ( JJJ.LT.I II-4)G0  TO  198 
III  1.2-11.1 12  *■  ABS  (  E 1  (J  J  J  )  ) 

19  3  CONTINUE 

ZZZZ1=EPS*(ABS(ZZZZ1)  +  1.) 

I F ( 11112. GT  .ZZZZI  > GO  TO  1 

ini  i=o." 

IIIL2=0. 

DO  199  JJJ=1,III1 
ZZZZ1=ZZZZ1  V'OfJjJ) 

I F ( JJJ.LT.I I1-41G0  TO  199 
11112  =  11112  «-  ABS  (  E2  (  J  J  J  )  ) 

199  CONTINUE 

ZZZZ1  =  EPS«(  4dS(  ZZZZI)  *■  1  .  ) 

I F ( 11112. GT. ZZZZI ) GO  TO  1 
ZZZZ1=0. 

117.12  =  0. 

DO  200  JJJ=1,  I  I  II 

Z 1 1 Li =1111 1  +  T  MP1( JJJ) . 

IF! JJJ.LT  .  I  I  1-4 )G0  TO  200 
l LI  12  =  III  12 _ +  ABS ( TNP1 ( JJ J )  ) 

200  '  "continue 

ZZZZl=EPS*( ABS(ZZZZl)  ♦  l.J 
IF( ZZZZ2.GT. ZZZZI ) GO  TO  1 
‘  ZZZZ 1=0. 

7.1112  =  0  . 

JO  20  1  JJ J=i, I  I  11 

'  ii in  =  ?. mi  Vtmp 2 (  jjj) 

IF (JJJ.LT  .1  I  1-4 >GC  TO  201 
11112=11112  *  ABS ( TMP2( JJJ )l 
u‘l  CG NT  IN JE 


1 


ZZZZ1=EPS 
I  F ( ZZZZ2. 
Z  ZZZ  1^0. 
Ill  12=0 . 
DO  202  JJ 
llllinil 
I  F (  JJJ.LT 
ZZZ  Z2  =  ZZZ 
CONTINUE 
i'll  11  =  EPS 
IF (  11112. 

_Gf  TO  2 
CUNT  I NUE  " 
CONTINUE 
00  203  JJ 


I  Ft  A  ti  S  <  K 1 
KKK= J J J 
GC  TC_2  04 
"CONTINUE 
Z ZZZ 1=0. 
KKK1=KKK 
DO  203  JJ 
11171=111 
I F (  ZZZZL/ 
DC  206  JJ 
I  F  (  ABSUl 
K  K  K  =  J  J  J 


GC  TO  207 

CONTINUE 

ZZZZ1=0. 


KKKl=KKK 
00  208  JJ 
ZZZZ1=ZZZ 
TFTzTZTl/ 

DO  209  JJ 
IFl ABSI91 


KKK=JJ J 
GG  TO  210 
CONTINUE 
ILL Z  1=0. 
KKK1=KKK 
DO^U  JJ 
'  Liiii =l:li 
I  F (  ZZZZ1/ 

n  C _ 2  12  _J  J 

I  F  (  A  G  S  (  K  2 
K  K  K  =  J  J  J 
GC  TO  213 
CONTINUE 
llll 1=0. 
KKK1=*KK 
UC  214  JJ 
L Llll  =  ILL 
I  F (  Llll  1/ 
'DC  215  JJ 
IF(AUS(~2 
K  K  <  =  J  J  J 
GC  TC  216 


MAOS(ZZZZl)  + 
GT.ZZZZDGO  TO 


J=l, I  III 

Z1  tTFAKEtJJJ)  _  _ 

.III “4 )G0  TO  202 
Z 2  +  ABS ( TFAKEt JJ J  )) 


MABS(ZZZZl)  ♦ 
GT.ZZZZ1 ) GO  TO 


J= 1*  Ml  II 


( JJJ) ) .LT.EPS)  GO  TO  203 


J=KKK1 ,NI I  I 
ZUA6SIRK  JJJ)  ) 

ABS( P 1 (KKK l ) . GE.  L ) lNO=IND 
J  = 1 , N II I 

(JJJ)  ) .LT.EPS)  GO  TO  2 Oo 


J=kKKl ,N  I  1  I 
Z  l+ABS^t  R_U  J  J  J  )  ) 

A  6$  (  R  1  (  K  KK  )  )‘7GE  .'Tf  IND=  1'ND 
J=1.  NI  II 

I JJJ) J.LT.EPS)  GO  TO  209 


J  =  KKK  1  ,  N  I  I  I__ 

Z 1  + Art's  CR  1I  JJJT7 

ABStR l(KKK) ) .GE 

J=  l,  NUI 

(J  J  J  )7  .LT.EPS) 


1 )  IND= I ND 


GO  TO  212 


J=  KKK 1, Mill 
Zl*AbS(R2(  JJJ)  ) 

ABS(  D-2(KKK)  )  .GE.  L  )  I  ND=  I  ND 
J= 1 , N I  1 1 

(JJJ) ) .LT.EPS)  GO  TO  215 


. *«**-***;*■  i  ' 


mm!' }mm 


215  CONTINUE 

216  ZZZZ1=0. 

_ KKK1=KKK  ♦  1  _ 

DO  217  JJJ=KKKl,NIlT 

217  ZZZZl  =  ZZZZl  +  AdS(R2I JJJ)  ) 

_ I F  (  ZZ  Z  Z  1  /  A  Q S  (  R 2  ( K KK  )  )_.  GE  .  1  )  I!MD=IN0  ♦  l 

OU  2  1 B  JJJ=1»NIII 

IFIA3SIR1I JJJ)). LT. EPS)  GO  TO  218 

_ KKK= JJJ _ _ _  _ 

GC  TC  219 . ” . .  ‘ . ~ . ~ 

213  CONTINUE 

219  _ ZZZZ1=Q. _ _ ___  _ _ 

K  KK  1  =  K  K  K  +  1 

OG  220  J J  J=KKK 1 »N 1 1  I 

220  ZZZZ1=ZZZZ1+ABS(R1IJJJ)  ) _ _ 

IFlZZZZl/A8S(Ri(KKK)).GE.l) INQMNO  +  1 
00  221  JJJ=1,NI II 

IFIABSITVSI JJJ}). LT.EPS)  GO  TO  221 


KKK* JJJ 
GC  TC  222 

III _ CONTINUE _ 

222  ZZZZ1=0. 

K  K  K  1 =  KKK  +  1 

_ DO  223  J  J  J  =  KKK 1 t  N I  I  I _ 

223  ZZZZ1=ZZZZ1+A3S(TVS( JJJ) ) 

IFI ZZZZ l/ABSITVSI KKK) ).G6. 1 ) I N0= IND  +  l 

_ DC  224  JJJ=1,NI  II _ 

IFI A 8S I R2 I J  JJ  ) ) .LT.EPS)  GO  TO  224 
KKK=JJ J 
GC  TC  225 

224  CONTINUE 

225  ZZZZ1=0. 

KKK1=KKK  +  1 

DO  226  JJJ=KKK1,NI  I  I"  ' 

226  ZZZZ1  =  ZZZZ1+ABSIP2( JJJ)  ) 

IF ( ZZ  ZZl/ABSI R2IKKK) ) .GE. 1 ) IND= I  NO  +  l 
DC  227  JJ J=l, NI II 

IF ( ANSI TWZI JJJ) ). LT.EPS)  GO  TO  227 
KKK- JJJ 
GO  TO  223 

227  CONTINUE 
223  ZZZZi=0. 

K  KK1  =  KK  K  +  1 

DO  229  JJ J  =  KKKl »N  I  I  I 

229  ZZZZl=ZZZZi+A8SITWZIJJJ)  ) 

IFI ZZ ZZl/ABSI TWZl KKK ) ).Gt. 1 ) 1 NO=INO  V" 1 
00  230  J J  J= 1 f  Nl II 

IFI  ABSU1I  JJJ)  )  .LT.EPS)  GO  TO  230 
KKK= JJ J 
GO  TO  231 

230  CONTINUE 

231  l  ILL  1=0. 

KKK 1 =KK  K  +  1 

DO  232  JJ J=KKK1 ,NI I  I 

232  ZZZZl  =  ZZZZl  +  ArlS(R  l(JJJ)  l 

IF  I  7ZZZ  1/  A.3SI  IU  (KKK))  .GE.  1  )  IND*  IND  *  1 
DC  233  JJJ=l»NI II 

1 F  f  A  3  5  I  3 1 1 J  J J ) ) . L T . CPS )  GO  TO  233 


KKK  =  J  J  J 
GO  TO  234 

CONTINUE _ _ 

Z ZZZ 1=0. 

KKK  1  =  KKK  4  1 

DO  233  JJ  J=KKKl,NI  I  I _ _ _  _ 

ZZZZ1=ZZZZ14ABS(R1< JJJ)  ) 

IF( ZZZZ1/  \BS(R1(KKK) J.GE.l ) IND=INO  +  1 

DC  236  JJJfU.ilI.IJ _ _ _ 

IF  (  ASS  f  R  1  (  JJJ)  i'.LT.EPS  I  GO  TO  236 
KKK=JJ J 

GO  _TU _ 2  3_7 _ _ _ _ 

CONTINUE 
till  1=0. 

KKK1=KKK  *  1 _ _ _ 

00  233  JJJ  =  KKKI»NI  I  I 
ZZZZ1=ZZZZ14A8S(R  l  (  J  J  J )  J 
I F  (  ZZZZ1/  A3S  (  R 1  (KKK  )  KGE_.  1J  IND=INO  4  I 
DC  239  JJJ=l,Nl  II 

IF( ABS( T V S ( JJJ) >. LT.EPS)  GO  TO  239 
KKK= J J J 
GO  TO  240 
CONTINUE 

ZZZZL  =  0« _ 

KKK 1  =  KKK  +  1 

00  241  JJJ=KKK1,NI  I  I 

ZZZZl=ZZZZL+AaS(TVS(JJJ)I 

I  F  l  ZZZ  Z  1/  AtJS  {  TVSI  KKK ) ) . GE. 1 ) I  ND=  I  NO  +  1 

DC  242  J  J  J=  1 « .N I  II 

IF ( A3S(R2( JJJ) ) .LT.EPS)  GO  TO  242 

KKK=J J J 

GC  TC  243 

CONTINUE 

ZZZZL  =0. 

KKK I =  KKK  +  1 

DU  244  JJ J=KKK1 ,N  1 1  I 

ZZZZl  =  ZZZZl  +  AdS  IR2UJJ)  J 

IF(ZZZZl/ABS(R2(KKK) ) .GE. I ) IN0= I ND  +  1 

DC  245  JJJ=1,NIII 

1  rl  A  BSTR  2  ( J  J  JTTTL  TTE  P  S )  G0~TCT"2‘4 5” 

KKK= JJ J 
GU  TO  246 
CONTINUE 
ZZZ Z  1  =  0. 

K  KK 1  =  KK  K  4  L 

OG  24  f~J  J  J*KKK  I » N TIT 

ZZZZl  =  ZZZZL  +  ABS(R2( JJJ)  ) 

IF(  Lilly/ A3SjR2(KKK)  )  .GE.  I )  IND=  I  NO  4  I 
DO  243  JJ J=1 ,Wi II 

IFC AHSIR2I JJJ) ).LT.EPS)  GU  TO  248 
KKK  =  J  J  J  _ _ 

GC  TO  249  '  ' . . . ~ 

CO NT  I NUE 
III  11  =  0. 

K  KK  1  =  K  KK  ♦  L 

DO  250  JJ J  =  KKKi ,N I  I  I 

ZZZZ1 =ZZZZI+4rtS(R2( JJJ) ) 

IP ( ZZZZ1/ABSI R2IKKK) ).GE.l) I ND=  INO  +  1 


00  251  JJJ=1,NI  II 

IF! ABSI TWZI JJJ) J.LT.EPS)  GO  TO  251 

_ K  K K  =  J  J  J _ 

GO  TO  252  "  ~ 

251  CONTINUE 

252  ZZZZl=0. _ _ 

KKK 1  =  KKK  +  1 

DO  253  JJ J  =  KKK1 ,NI  I  I 

_25  3 _ l  III  1  =  ZZ  III  *  A  B  S IT  U  Z  I  J  JJ  )  ) _ _ 

I  F ( ZZ Z Z 1 / AbS I T mZ( KKK )  ) • GE . 1  )  I N0= I NO  ♦  1 
NI  I  I  =  NI  I  I  *■  1 

_ ZZZ31  1)  =R  II  NI  1 J  J _ _ _ 

DO  254JJJ=2»NII  I 

254  ZZZBI  JJJ)=KHNIII-JJJ«-1)  +  T*  ZZ  Z  B  (  J  J  J-l ) 

YNEWI  1 ) =Z  ZZ  B ( NI 1 1  ) _ 

Z Z Z 3 ( 1) =R2 I NI I  I  ) 

DG  255JJJ=2,NI 1 1 

255  ZZZ3( JJJ ) =R2( NI II -JJJ+1 )  +  T* ZZ Z B C J J J-l ) 

YNE'wl  2)=ZZZB(NI  II  )  '  ‘  . . 

ZZZB( 1 ) =XH1 ( NI  I n 

_ DG  2  56JJJ=2fNI  I  I _ 

256  ZZZBI JJJ)=XMIINII I-JJJ+1)  +  T*ZZZB( j J J-l) 
YNEWI  3 ) =Z  ZZ  B( NI II ) 

_ Z  ZZ  3  (  l)=XM2(Nim _ _ 

DO  257J JJ=2,NI I  I 

257  ZZZBI JJJ)=XM2(NII I-JJJ+1)  +  T*ZZZ3IJJJ— 1) 
YNEv»(  4)  =ZZZB(  NI  II  ) 

ZZZBI  1) =  E 1 ( NI II ) 

DO  258JJJ  =  2*NI I  I 

253  ZZZBI JJJ)=E1(NI II-JJJ+1)  +  T* ZZZBI J JJ-1 ) 
YNEWI 5)=ZZZB(NIII ) 

ZZZBi  1 ) =E2 I NI  I  I  I 
DC  2  59 J  J J  =  2  »  N  1 1 1 

'259  ZZZBI JJJ )=E2(NI  1 1  -  J  J  J  + 1 )  +  T*ZZZ3( J J J-l ) 

YNEWI 6)=ZZZ3INIII I 
ZZZ3I 1 I =TMP1( NI II ) 

DC  260 J J J =2 » N I  I  I 

260  ZZZBI JJJ) =  TVP1 INI II-JJJ+l)  +  T^ZZZBI J JJ-l ) 
YNEWI 7)=ZZZBINI II ) 


ZZZBI  U=TYP2(  NI  II  ) 

DO  261JJJ  =  2,NI  I  I 

261  ZZZBI JJJ)=TMP2(NI II-JJJ+1)  +  T*ZZZB( J J J-l ) 
YNEWI  3)=Z  ZZdlNl  II  ) 

ZZZBI  l ) =  TFAKE (Nil  I ) 

DG  2o2JJJ=2tNI  1 1 

'262  ZZZBI  J  J  j )  =T  FAKE  IN  1 1 1  -J J J+l )  +  T*ZZZB(  JJ J-ll 
YNEWI 9 ) =ZZZB( NI II ) 

RE  TJRN 


APPENDIX  B 


Output  of  SERIES  for  constant  spacial  derivative  formulation 


! 

\ 


SUBROUTINE  C  I  f  FUN  (  T  » Y  »0Y  ) 
DIMENSION  0Y<20),Y(20I 

fl  A  T  A  0/1 


D Y ( I ) -0 • 

DY(2) =-Y( 5) *Y( l2)-( YI5J-YC4 ) ) /JX 


0Y<4) =0. 

CY(5)=-Y( 5)*Y(5)/Y(2)*Y( 12)-<R*Y< lO)-<  Y(5)/Y(2) )»*2)*< YI 2 )-Y<  II)/ 
Y(2)*(Y(5)-Y(4))/0X  _  _ 


DYt  6)*-YI 6l*Y(o)/Y(3)*Y( 13 J - ( P»Y ( 1 1 1- C Y ( 6 ) /Y ( 3 ) J **2 )* I Y ( 3 )- Y ( 2 J ) / 
*DX-2.*YI6)/Yt3I*lYI6)-YI5)I/DX 
HYt  71=0.  _ _ _ _ 


JY(8)*-Y<  5  I’M  Y(8I  /  YI2I+R*  Y(  10)  )*;  Y(  12) +Y (5 1 *Y( 8  I /Y<  2 1**2* (  Y(2)-Y(  1 
+  )  )/DX-(Y(e)/Y(2)+R*Y(10))*(Y(5)-Y(4))/DX-Y( 5)/YI2)*(Y(8)-Y(7))/DX 
Y(6)MY<9)  /  Y  (  3  )  +  R  *  Y  (  11  ) )  *  Y(  13  ]±  Y  (  6  )  »  Y  (  9  )  /  Y  t  3  )  »*2»  (  Y  (  3 )  -  Y  t  2 


♦  )  )/DX-<Y<  9) /Y  (3  H-R*Y  { 11  ))*(  Y<  6)-Y  (5)  l/DX-Y  (6  >/Y(3»  MY  (9  )-Y(  8)  )/CX 
DY<  10=0. 

D  Y  (  1 1 )  =  0 


CY ( 12 )=0 
DY< 13) =0 
RETURN 


ENC 


SCaRCUUNE  SCUT,  YCtYNEWt  INC) 

OIPENSICN  YC(20),YNEW(20) , ZZZBI 20 ) , OROl 20 ) ,RO( 20) , DR  LI  20) ,Rl(20> , 
♦XK1(20),XK1(20).TST(20)  ,  TSC  <  20) ,XMO  (.20 )  ,  J  SV  (  20)  ,TSW(20J  ,TSX(20I  .0 
♦R2I20) ,R2(20) , XP2 I  20) , XK2 1  20 ) , TSY 1 2 0) fTSZ(20l ,TS0( 20) #TSi (20) ,TS2 
♦(20) ,DXM0(20)  ,CXMU20) ,TS3( 20) ,T$4( 20) ,TS5(  20) .  TS6 1 20) , TMP 1 1  20 ) , T 
♦S7(20).TSE(20).7S9t20).ITS(20).TTT( 20) , T TU( 20 ) » T T V < 20) . TTw ( 20) . TT 
+  X(20)  ,  TTY  120  , TT0I20) ,  TTK2C)  ,TT2(  20)  ,DXM2(  201 ,TT 3 1  20) , T T4 1  20 ) , T T 
♦ 5( 20) ,  TT6 (20) , TMP2( 20) ,TT 7(20) ,TT8(  20) ,TT9(20) ,TUS(20) ,  TUT  t 20) ,TU 
♦  01  20)  jTUV(20)  .TOW  (  20)  ,  TUX  120  ,TUY(  20)  ,TUO(  20)  .  TUI  120) ,  TU2  1 20)  ,  DEO 


♦  (20) , E0( 2  G ) » OE 1 (20 ) *  E 1 ( 20 ) »TU3(20) , Tub ( 20 ) »  TU6 ( 20 ) »TU7( 20) »  TUfi 

♦  (  20)  ,  TU9(20  )  ,TVS(  20)  ,TW(20)  ,TVVI20)  .TVWI20)  ,TVXI  20)  ,TV2<  20)  ,TV3 

♦  (  20)  .  TV4(  2  0  ) »  T  V6  (  20  )  ,T  V7(2C)  ,  TV  8J2  0  >jj  V9 1 2  0)  .DE2(20)  .62  (20)  .TwS 
♦(20), TWC( 20 ) , TWV( 20) ,TWW(20> ,TWXl20) , TWYI 20), TWZI 20) ,  TWO! 20) ,Tw2 

♦  (20),Tw3(20),TI»4l20),TW9(20)»TXS(20),TXT(20)»TXV(20)tTXw{2GI»TXX 

♦  (20)  «  TXV  (  20)  ,  JXKK20)  ,  DXK2  (  20  ) ,  D  TMP 1 ( 20 ) , DTMP 2 ( 20  ) _ 

DATA  R/1716./,CX/.01/ 

EPS=1.0£-6 
R0(  I )  =YC  (  1) 


R 1 ( l ) = YC( 2 ) 
K2 (  L)=YC (3) 
XI-OI  L  ) =  YO  (4  ) 


XP1(1)=Y0(5) 
X P2 ( 1 ) = YO (6  ) 
EO ( 1 )  =  YO  ( 7 ) 


E 1  (  1 ) =  YG  <  €  ) 

E2(  1)=Y0(  S) 
TPPll 1)=YC(  10) 


TMP2I 1 1  =  VC (  11) 
XK1 1 1 ) =Y0 ( 12 ) 
XK2( I )=YC( 13) 


IAC  =  0 

OC  1  111=1,19 
N I  I  I  =  II  l 


1111=111  -  1 

IF  t  I1I.GT.  I  )  C  P-  0  (  I  I  I  )  =0. 

QRC( 1 ) =0, 


R 0 ( 111  +  1 ) =QRO( I II)/ FLOAT!  Ill) 
T  S  T ( 1 II ) =0  . 

OC  100  JJ J=1  ,  I  I  I 


T ST (  l  1 1 )  *TST  (  1 II)  ♦  XML  I  JJJ)*XK1(  1 1  l-JJJ*l) 
rsu( III ) *«T ST ( 1 1 1  ) 

tsv(  i  m  =  xpi(  i  n)-xMOi  mi _ 

TSW(I1I)=TSV(  IID/OX 
T  £  X  {  1 1  I  ) =TSO<  l  m-TSW(  III) 

DPI ( I  1 1 ) =TSX(  1  1 1) 


R 1 (  ill  ♦  1)=0R1(III)/FLCAT(  III) 

T  SY ( I  II ) =0. 

DC  101  JJJ=1. 1 1  I  _  _ 


T  SY(  1  1 1 ) =TSY( 1 1 1)*XM2( JJJ)«XK2( I  II-JJJ*1) 

TSZ( I  II )=-TSY(  III  ) 

TSOI  I  I  I  )=XP2(  1  I  1  l-XMK  I  II)  _  _ 


tsk  iii»  =  rso(  im/cx 

F  52  (  III)  =  I5Z(Ili)-T£imi) 
0R2( 1 1 l)  =  fS2(  i  1 1) 


TS3Uin=C. 

DC  102  JJJ=l,lII 

T  S3  (  I  1  1 ) =  T  S3  (  I  1  I) ♦XMII JJJ)«XP1(  1 II-JJJ+  1 ) 
IF  (III. EC-1)  GC  70  103 
T$4(  III)  =  TS3(  I II)-TS4( l)*Rl(l  II) 

IF ( I  I  I  .EC. 2  )GC  TU  104 _ 

DC  105  JJJ  —  2V 1 1 II 

105  TS*<  I  1 1  )=TS4<I  II)-TS4C  J  JJ)*R1(  II  I-JJJ+i) 

104  TS 4  (  III  )  =  T  S4  (  I  1 1 )  /fill  1) _ 

GC  TC  106 

103  TS4( III )=TS3( II I)/Rl( III) 

106  CONTINUE _ 

T  S  5 ( 1(11  =  0. 

DC  107  JJJ=  l,III 

107  TS5(im  =  TS5<iii)  +  rs4(jjj)*xKnm-jjj*n 
T$6< I  1 1 ) =-T S5 (  I  II  ) 

T  S  7 ( 1 1 I )  =  TPP1 ( 1XI)*R 

_ IF  ( i  II.EC.l)  GC  TC  108 _ 

rsat  in  )=xpk  1 1 1 ) - TS8 (  i>*rii  in  ) 

I F (  1 1  I . EC  .2  »GC  TO  109 

_ l)C  110  J  J  J  =  2  f  1  I  11 _ 

110  TSd(  III  )=TSE(  I  II)-TS8<  JJJ)*Ri(  1 1  l-JJJ«-l) 

109  t  sau  1 1  »  =  Tseu  i  n  /rki  » 

_ GC  TC  111 _ 

10a  T  sa  (  I  1  I  )=XIU  (  1  1 1)  /R1  (  1 1 1  ) 

111  CONTINUE 

_ TSS( III )=C. _ 

OC  112  JJJ  =  1, I  I  I 

112  TS5(  III  )=TS9(  III)+TS8(  JJJ)*=TS8(  I  Il-JJJ*l) 
TTS( I  1 1 )  =  TS 7 ( 1 1 1 » - T S9 (  I  I  I  I 

TT7( III)=H1 (lli»-«0( III) 

T  TC ( 1 1 1 ) =0  . 

_ OC  113  JJJ=1,  I  I  I _ 

113  TT'Jl  I  1 1  )  -  T  T  L  ( IlII+TTS(Jjj)*TTT(  I  Il-jjj+l) 
TTV( 1 1 1 ) -TTU( 1 II) /CX 

T  T  (  1 1 1 )  =  T  Sfe  1  1 1  I  )-TTVl  I  11) 

T  T  X (  1  1 1 ) - XP  1  (  I  1 1 ) *2. 

IF  ( 1 1 1 .EC.l)  GC  TC  114 

_ TTY(nn«TTX(llI)«TTYIH»RHHl) _ 

IF( II I.EC.2IGC  TO  115 
DC  116  JJJ=2t I  III 

116  TTY(IH)  =  7TV(111)-TTY(JJJ)»KKII1-JJ  J <- 1 ) 
115  T  T V( I  1 1 )  =  TTV ( 1 1 1) /Rl(l) 

GC  TC  117 

114  TTV(  I  I  I  )=  7TX{  1  ID  /  R  l  (  Ill) _ 

117  CCMINOF  . . 

T  TO ( I  I  I ) =0. 

_ uc  tie  jjj  =  i ,in _ _ 

113  TTO(  1  I  I  l*T  T0~(  II II  *TTY(  JJJ)*TSV(  l  ll-JJJm' 
Til ( I i I ) ®TTOC 1 1  I) /CX 

_ IllLLLLLi-LT M  III)  —  TTI(III) _ 

uXR  l(  I  III  =t"72  (1  II) 

X/llIII  ♦  l)  =  CXRl(llI)/FLCM(lll) 

_ TT3(  1  m  =C. _ _  _ 

OC  119  JJ J=1 ,  I  I  I 

119  T  T7I  1  1 1  )  =  TT3( 1 I II+XM21 JJJ)*XH2( I  1 I-JJJ+l) 

_ IF  (HI  .EC.  1)  GC  TC  1 2 0 _ _ 

Tl4( I  I  I  I  =  T T 3 (  I  1 I )  •  T  T4  ( l)*M2<  I  II ) 


IF l III .EG.2 )GC  7 U  121 
OC  122  J J J  =  2* 11 II 

122  TT4(1I1)=TT4(II1)-TT4(JJJ)»R2(I1I- JJJ»1  ) 
121  TT4<  III)=1I4<  I1D/R2U) 

GC  TL  123 

120  TT4(11I)=TT3<1I1)/R2(III) _ 

123  CONTINUE 
TT5Uin*0. 

_ CC  124  JJJ=1,1I1 _ 

124  TT5(  1 1  l ) =T  T  5  l  1  i  I)  *TT4l  JJJ)*XK2(  1 II-JJJD) 
7  T6(  I  ll»  *-TT5(  i  ID 

_ TT7  I  I  1  I  )  =  T  MP2 (  1  11  HR _ 

IE  I  III .EC. 1)  GC  IC  12  5 

TT3< II I»  =  XM2( 1 1 I)-TT8< l)*R2(  I  III 

_ IF(  III.EG.21GC  TO  126 _ 

DC  127  J J J  =  2  »  I  1 11 

12  7  TTe(III)=TT8( 1 1 1) -TT8 I J J J ) *R2 ( 1 1 1- J J J* 1 1 

120  T  T8 ( 1 1 1 ) =TT8<  I1I)/R2<1) _ 

GC  TC  123 

125  TT8I  11U-XH2I  11I1/R2I  III) 

123  CONTINUE _ 

T  T9 ( I II1=C. 

DC  129  JJJ=l,III 

129  TT5( I  II )*TT9( I II) »TT8( JJJ)«IT8( I  I I-JJJ+1) 
TCS(III)-TT7{II1)-TT9(  III) 

TUT(  III  )*R2(IIl)«Rl(Uil 

_ TliU  II1II-C. _ 

DC  130  JJJalflll 

130  TULUII)*TUUf  III)*TiJSUJJ)*TUT(IlI-JJJ*l) 

_ TUV  (III)  gTUU(  im/cx _ 

T  UM  1  1 1 1  *TT6  (  III)*TvJV(  III) 

TUX  ( 1 1 1  ) =XP2 (  1 1 1 )  * 2  . 

_ IE  (  III.EC.1J  GC  TC  131 _ _ 

TUY (  1  II)  =  I UX { l l 1 1  - TUY  C 1)  *R2<  I  II  ) 

I F  (  1 1 l .EG .2  )GC  TO  132 
_ OC  1 3 J  J JJ  =  2 1  I  III _ 

133  TUY{  1 1 1  )  =  IUYI  l  i  l  )«TUYUJJ)*R2<  I  l  l-j'jJU  I 

132  TUY( III)*7UY( I  II) /R2(i) 

_ GC  TC  134 _ 

131  TJYl I  1 1 )  =  TU  X ( II  I ) / R2 (  I  II) 

134  CONTINUE 

_ TUQ(  1 1  IDO. _ _ 

DC  135  JJJ=1,  I  I  I 

135  TUO( I  1 1 ) =TUO< 1 1 U +  TUY<  JJJ) *TSO{  I  I I-JJJ+1) 

_ TUI  (  1  I  I )  =  TUO C  II  IJ  /  OX _ 

TU2 (  l II ) =TUw(  lil)-TUl(III) 

OX V2 ( 1 1 1 ) =  TC2 (  1 1 1  ) 

_ X  y  2  (  III  +  1DUXH2I  I  m/FLCATI  II  IJ _ 

IF  (  III.GT.l ) 0£C(  III  J=C. 

OtO(  1 ) =0 . 

_ EOj  II  I  *  1) =  DEO < I  I  I  )/FLCAT (III) _ 

IF  (  I  II.EC.~1)  GC  TC  136 
TJ3III  I)=E1(  1 1 1  )-TU3U)  “R1I  II  I) 

_ 1  F ( 1 1  I  . EG  .2  ) G G  TO  137 _ 

JL  13  8  JJJ  =  2,  II  II 

133  T:J3(I  ID  =  TU3<  1II)-TU3( 

137  Tui(IIl)  =  TL3(I  U ) / K  1  (  1 ) _ 

GC  TC  139 


1J6  Tl>3(  1 1 1 ) =  El  (  1  2II/R1CII1I 

139  CONTINUE 

TC5(11I)-TU3(  HI)»TS7(1I1) _ 

TU6 ( l I  I  1 =0  . 

0C  190  JJJ  =  1 ,  1 1  1 

140  TUt»(  II  1  )  «  TU6(  im+Xi'U(JJJ)»lU5UH-JJJ»l  ) 
TJ7(  im=c. 

DC  141  JJJ  =  1,  I  1  1 

191  T  t  J  7  <  H  I  )  =  T  Li  7  (  1  II)  »  T  U6  (  J IJJ)»XK1(  1  I  I-JJJ*  U 

Ti-aiiii  » =  - t u 7 <  i  ii ) 

TU9U  m=c. 

PC  LjuL  JJJ=l»  HI _ 

1*2  TU9U  II)  =TL9(  I  1 1)  *XMK  JJJUEU  II  I-JJJ+1) 
TVS!  I  I  n=o. 

_ PC  143  JJJ=l.  II  I _ 

143  TVSU  l  l)  =  T  VSI  I  I  l)  *fU(  JJJUUlI  1 1  I-JJ  JU  ) 

IF  (  III. EC. 1)  GC  TC  144 

_ TVTI1H)=TU9(1II)-TVT(1)«TVSUII) _ 

IF l  1 1  I. EG. 2  )GC  TO  14b 
DC  146  J JJ  =  2  »  I  1  1 1 

196  T  V  T  (  III)=TVT(  III)-TVT(JJJ)»TVS[  III-JJJ+1) 
145  T VT ( 1  1 1 )  =  TVT< Ill) /TVS(l) 

GC  TC  147 

144  TVT(  I  I  1  )  =  TU9(  1 1  1) /TVS!  I  III _ 

147  LCMINUE 

T V V (  I  1 1 ) =0 . 

_ DC  148  JJJ  =  1 , 1  I  I _ 

148  TVVI  III)=TVVUli)*TVT<  JJJ)*TTT(  III-JJJU) 


Tv*<  III)  =TVV  C  IIU/DX 

_ ivxiim»Tuec  nn+ivw <  1 1 n _ 

rv2i ii i ) =o. 

DC  149  JJJ=l,  I  I  I 

149  TV2(  1  II  )  —  T  V  2  (  II1)+TU5(  JJJ)MSV(  II1-JJJ+1) 
TV3(  I  I  I )  =  T V 2  t  I l I > /OX 

TV4 ( I  I  I ) =TVX(  I  I  I ) -TV  3 (  III) 

_ T  V 6 (  I  I  I  i  =E  1  (  I  I  I  )— 601  HI) _ 

TV  7 ( I  I  I ) =0. 

DC  150  JJ J=l,  1 1  I 

150  Tv7(nn  =  Tv7<nn»TS8(jjj)«TV6iin-jjj»i) 
1V8I I II)  =  TV7( I  1 1) /OX 

TV9I  I  1 I )=TV4<  1ID-TV8U  1 1  ) 

_ DE  1  I  I  I  I  )  =TV9(  III) _ _ 

E  1  U  l  I  ♦  l)*l)El  (I  1 1  l/FLOAT  I  1 1  I) 

IF  ( 1 1 1 .EC. 1 )  GC  TC  151 

_ TwS (  III )=E2 II1I)-T^S(I)*R2(III) _ 

l F  I  i l I . t w .2 ) GC  IQ  152 
PL  153  JJJ*2»  II  II 

153  TVo!  1 II  )  =  TaS(  I  H)-TWS(  JJJ  I  >K2(  I  1  I-J  JJ  ♦  IJ 
152  Tv.SU  I  I)  =  TW5m  D/R2I  l") 

GC  TC  154 

151  TUS( I  I  I )  =  E2( 1 1  1  )/R2(  1 1  1  ) _ 

194  CCNTIiluE 

TrtUl  1 1  I  )  =  T  V*  S  (  IID+TT7I  III) 

_ TwV(  I  1 1 ) =  C  . _ 

DC  15  5  J  J  J  =  1 ,  I  i  1 

13J  Tv.  VI  I  1 1  )  =  r*V<  III )  ♦X.YZI  JJJ>*T»U<  I  I  l- JJJ* II 

_ T  «  4  (  I  I  I  )  g  C  . _ 

uC  156  JJ J=  l ,  i  I  I 


15o  Ti*w(lil)-TV'W(III)*TMVtJJJ)*XK2(III~JJJ  +  l) 

TkXnil|>-ThW(lII) 

_ T  h  Y  I 11II»0. _ 

DC  157  JJJ=1,  III 

157  TWYI  1  l  I  )=Tl*Y(  I  1 1)  +XM2  (  JJJ)*E2(  1 1  l-JJJ+U 

_ n»/.<  i  m  =0. _ 

DC  158  JJ J=  l,III 

158  Ti*Z(IIU*TkZUlIl  +  R2I  JJJ  I  +  P21  II  I-JJJ+1) 

_ IF  (  I  ll.EC.l)  Gl  TC  1 55 _ 

TinO  (  I  II  I  =  Tl*  Y  <  I  I  l)-ThO(iJ*Tw2(  111) 

IF( 11I.EC.2IGC  TC  160 

_ JC  lfcl  JJJ  = 2,1111 _ 

161  T«0< 1 1  I )=T*0< I  I l)-TrtOI JJJ)*TW2( I  II-JJJ+1) 

160  T In 0 (  lll)=TWO(III)/TwZ(l) 

_ GC  TC  162 _ 

15‘i  TlrtO  (  11I)  =  TWY(  1 1  I)  /  T  wZ  (  Ill) 

162  CONTINUE 

_ Ta2( 1 1 1 ) =c « _ 

CC  163  JJJ=l,  III 

163  Tk2(  1 1 1  )  =TW2(  1 1 1 )  *TWO(  J  J  J  J  -« TUT  C  I  II- JJ  J+i  1 

_ T»3( III ) =TU2 ( III) /CX _ 

Th4U  1 1)  =TWX(  II I)  +TW3I  HI ) 

TW5( l l l )=0. 

_ DC  164  JJ J  =  1  ,  I  n _ 

164  TWSII  UI  =  TW9(U  II+TWUI JJJ)»TSOl  I  II-JJJ  +  U 
TXSI  I  II  l«T'«S(  IID/OX 

_ TXT!  HI)  -T>>41  11  n-TXS<  III) _ 

TXVl  I  I  l)=62U  m-El(il  M 
TX a  I  I  1 1  )  =0  • 

_ DC  165  JJJ=1,  111 _ 

165  rxwf I i I )  =  IXW| 1 1 1) ♦TT8(JJJ)*TXVl  I  II-JJJ+1) 
TXXI  I  1  I ) =  TXw< 1 1  1) /CX 

_ txyi  i  i  ii  =  txti  i  m-Txxnin _ 

DE2(  I  I  I) =TXY( II  I) 

E2IIII  «■  1)=UE2(I I l)/FLCAT<  III) 

_ IF  (  lII.GT.l)UTfrPl(HI)=Q. _ 

OTF<Pi<l)=C. 

TKPK1II  ♦  l)  =CTMPl(  II  M/FLCATIl  II) 

_ IF  (  III.GT.IIDTIYPZI  IIIHJ. _ 

DTPF2 ( 1 ) =C« 

T iM P 2 (  ill  +  1 1 *U IMP2 (  I  III/FLCATI  III) 

_ IF  (  I  1 1 ,GT  .  1  )CXK1 (  III ) =0. _ 

DXK 1 ( 11=0. 

XKl(  111  +  U=UXKi(  lin/Fl.CAT(  III) 

_ IF  Illl.GT,  1  )DXK2  (  II  H=0. _ 

0XK2I l) =0. 

X K 2 (  III  l»  =  CXK2(im/FLCAI(  III) 

_ IF  ( 1I1.LT  .<*)GC  TU  1 _ 

1111=111  ♦  l 
ZZZZ1=0. 

_ ZZZZ2=0« _ 

DC  16b  JJ J  =  l,  I  I  11 
ZZZZ1  =  ZZZZ1  0 ( J  J  J  ) 

_ 1  F( JJ J.LT  .1  1  1-4 )Gl  TO  166 _ 

Lt  L12.-LLLU  *  ABSIPOl  JJJ)  ) 
l6o  CONTINUE 

_ ZZ  LL  1  =  CPS*  (  AcS  (  ZZ  Z Z  1 )  ♦  1 .  ) _ 

IN  ZZZZ2. GT  .ZZZZIJGU  IL  1 


_ nc  16  /  jj j=  1, 1 1  n _ 

,  '  ZZZZ1=ZZZZ1  ♦RllJJJl 

'  IF( JJJ.LT.I II-41GC  TG  167 

_ LLILI  -  LLLL'l  »  AdS  <R1( J JJ)  I 

16  7  CCMINUt 

^  ZZZZ1=EPSM AdS(ZZZZl)  +  1.) 

_ IK  ZZZZ2.GT  .ZZZZ1  )GC  TC  1 

ZZZZ1=0. 

Z  ZZ  Z  2  =0 . 

_ OC_  16 d  J J J=  It  1  Ill _ 

ZZZZI=ZZZZ1  *R2 ( J  J J  1 
IFIJJJ.LT.  I II-41GC  TO  168 

_ ZZZZ2  =  ZZZZ2  ♦  A  8S  (  R2  ( J  J  J  1  1 

ltd  CONTINUE 

ZZZZ1=£P$*(A8S(ZZZZI)  ♦  1.1 

_ IFIZZZZ2.GT.ZZZZ11GC  TC  1 

.  ZZZZI=0. 

ZZZZ2= 0. 

PC  169  JJJ=1,  1  1  II _ 

ZZZZ1=ZZZZ1  +  XPGI  J  J  J  ) 

IF< jJJ.LT.I II-41G0  TG  169 

_ ZZZZ2=ZZZZ2  +  ABS( XPOI J JJ ) ) 

lo9  CONTINUE 

ZZZZl=EPSM  ABSIZZZZl)  ♦  1.1 

_ l  h  l  ZZZZ2.GT.ZZZZ1 1GG  TG  1 

ZZZZ1=0. 

ZZZZ2=0. 

_ oc  HQ  jjj  =  ik  LU _ 

ZZZZ1=ZZZZ1  ♦XMKJJJ) 

IF ( JJJ.LT.I  II-41GC  TC  170 
_ ZZZZ2=ZZZZ2  *  AEStXMK  JJJj  1 

170  CONTINUE 
ZZZZl=EPS*(AbS(ZZZZl)  ♦  l.) 

_ IF  1 ZZZZ 2.GT .ZZZZ1 ) GO  TC  1 

j  ZZZZ1=0. 

ZZZZ2=0. 

f  _ CO  171  JJJ=1,  Till _ 

ZZZZKZZZZ1  *XP2(JJJ) 

IF ( JJJ.LT.I 11-4 1G0  TC  171 
_ UIL2  =  LLLll  *  AeS(XP2(  JJJ1  ) 

171  CONTINUE 
ZZZZl=EPSM/»eSIZZZZl>  ♦  1.1 

_ 1FIZZZZ2.GT.ZZZZ1 IGU  T C  1 

ZZZZ1=0. 

ZZZZ 2=0. 

_ 00  1  72  JJJ=  1  Kill _ 

ZZZZ1=ZZ.ZZ1  *  EO  (  J  J  J  ) 

IF( -JJ.LT.l  I  1-4 1GC  TO  172 


ZZZZ2  =  ZZZZ2  ♦  AOSIEOUJJU 


IMJJJ.LT. 1  11-4  )GC  TO  173 
11112=11112  *  ABSIEIIJJJ)) 

173  XCNTINU£ _ 

ZZZZ1=EPS*< A6S1ZZZZ1)  ♦  l.) 

IF  (  ZZZZ2.GT.ZZZZUGC  TC  1 

_ ZZZZ1=Q. _ 

11112=0. 

OC  174  JJJ=l,IIIl 

_ ZZZZ1-ZZZZ1  ♦  E2 ( J  J J  1 _ 

IF < JJJ.LI .1  1  1-4 1GC  TO  174 
III  12  =  l L 112  ♦  /» eS  (  E2  (  JJJ)  ) 

174  CCM1NCF  _ _ 

ZZll  1=EP$*  (  AOSt  ZZ<.Zl »  4  l.) 
IFIZZZZ2.GT.ZZZZL1GC  Ic  1 

_ ZZZZ  1  =  0. _ 

ZZZZ2=0. 

OC  175  JJJ=  1,  I  I  1 1 

_ lllll-lllll  T  N  P  l  (JJJ  I _ 

IFIjJJ.LT.  11  l  “4  )GC  ro  173 
III  12  =  111.12.  *  AES ( TPPI < JJJ )> 

173  CCM1NLE _ 

ZZZZi  =  EPSM  AES(ZZZZl)  ♦  1.) 
1F(ZZZZ2.GT.ZZZZ1  )GC  TC  1 

_ iim  =o. _ 

11122=0. 

DC  176  JJJ=1,1I11 

_ ZZZZUZZZZ1  +1PP2<JJJ) 

IFlviJJ.LT.I  1I-4JGC  TO  170 
11112  =  11112  ♦  ABS( TyP2( JJJ) ) 

176  CONTINUE _ _ _ _ 

ZZZZl  =  EPSMAdS~<ZZZZl)  +  1.) 

IF( ZZZZ2.GT .ZZZZl 1GG  TC  1 

_ ZZZZ1  =  Q. _ _ _ 

11112=0. 

DC  17  7  JJ J=  1, 1 1  II 

_ iiin=inn  »xk i ( jjj j _ 

IF  t  JJJ.LT  .111-4  »GC  TC  177 
11112  =  11112  *  AbS( XK1( JJJ  )  ) 

177  CONTINUE _ 

ZZZZl=tPS* ( ABS(ZZZZl)  ♦  1.) 
IF1ZZZZ2.GT.ZZZZ11GC  TC  1 

_ llll 1=0. _ 

11112= 0. 

CC  l 78  JJJ=  1,  I  III 

_ llll  1  =  ZZZZ  1  + XK  2JJJJJ _ 

IF(  JJJ.LT.  1  11-4  )GC  T C  17~8 
7.1112  =  12.112  *  AES<XK2< JJJ ) > 

178  CONTINUE _ 

ZZZZl=cPS* ( AdS( ZZZZl)  ♦  l.l 

I F ( ZZZZ2.GT. ZZZZl 1G0  TC  l 
_ CL  TC  2 _ 

1  cLM  INUF 

2  COM  I  MJE 

_ DC  17  7  JJJ=  UMI I _ 

IF l AuSIKl (JJJ) )  .L  T  .  EPS  1  GC  TC  179 
K  K  K  =  J  J  J 

_ GC  TC  180 _ _ _ 

1  7  i  CCM  1  DUE 


ISO  ZZZZ 1=0. 

KKK  l  =  KKK  +  1 

_ uc-  18  L  JJJ  =  KKK1,NH1 _ 

LSI  ZZZZ1=ZZZZL*AUS(R1 { JJJ)) 

IF  (  ZZZZ1/A3SI  Kl(KKK)  )  .  G£  .  1  )  INC  =  l  NO  *■  1 

_ nc  162  JJJ= 1  ,.NI_Ll _ 

If-  lABS(RKJJJ)  J.LT.EPS)  GC  TC  L82 
KKK  =  J  J  J 

_ GC  1C  18  1 _ 

132  CCM1NUE 

193  ZZZZ1=0. 

KKK 1=KKK  j  1 _ 

DC  184  JJ J^KKKl  ,NI  1  1 
13-*  ZZZZl  =  ZZZZl  +  AdS(Rl(  JJJ)  ) 

_ 1  F  (  ZZZZl/AbSIKltKKK)  )  .  G  C  .  I  )  I N  C  =  I ND  *  1 

DC  195  JJJ=l,Ni  (1 

l F t ABSIRltJJJ)  I.LT.EPS1  Gu  TO  185 

_ KKK=J J J _ 

GC  TC  166 

185  CONTINUE 

186  ZZZZ1=0. _ 

KKK 1=KKK  +  1  - . 

DC  137  JJJ  =  KKK1  »NI  I  I 

137  ZZZZ1=ZZZZ1+A J S ( R  1  ( J J J  )  ) _ 

IF<  ZZZZ1/ JHS(rT(KKK  H.CE.l  1  INDUING  «•  1 
DC  198  J J J=  1 1  M  11 

_ IK AES(R2( JJJ)) .IT. EPS)  GC  TC  188 

K  K  K  =  J  J  J 
GC  TC  1 8*5 

133  CONTINUE _ 

ld9  ZZZZ1=Q. 

K  K  K 1  =  KKK  +  1 

_ DC  190  JJ J=KKKl  ,N  1  I  1 _ _ _ 

190  ZZZZ1=ZZZZUa'6S(R2<  jjjm 

I  F  (  ZZZZ1/  4QS(R2(KKKJ  )  .  GE  .  1  )  IN.  C=INO  ♦  1 

_ DC  191  JJJ=1>N11I _ 

IF(AHS(K2(JJJ)")  .LT.EPS)  GC  TC  191 
KKK= JJ J 

_ GC  TC  192 _ 

191  CONTINUE 

192  ZZZZ1=0. 

_ KKK1  =  KKK  *  1 _ _ _ _ 

DC  153  JJ J  =  KKK1 ,N  I  I  I 
19  3  ZZZZI  =  ZZZZ1+APS(R2(JJJ) ) 

_ i f ( zzzzi/3 esjR2 <  k k km.ge .in m:=  ind  «•  1 

DC  19 A  J  J  J =  1 , N  I  I  I 
I  FI  ASS! K2 (JJJ) ) .LT.EPS)  Gu  Tu  194 

_ KKK  =  JJ  J _ 

GC  TC  195 

194  CONTINUE 

195  ZZZZ1=0. 


KKKl=KKK  +  1 

OC  159  JJ J  =  kKK1  ,N  I  i  I 


199  LLLLLiAAlJJ-1* jisou  (  jjj)  ) _ 

IF  (  ZZZZl/Ai3$(R  I  (  K  K  K  )  )  .GE.l  1  I N  C  =  1  N  Q  ♦  l 
OC  200  JJJ— If  i\  111 

I  F  1  Mii  TVS  t  J  J  J  )  ) .  LT  .EPS  )  GC  TC  200 _ 

K  K  K  =  J  J  J 
GC  TC  201 

CCMIMlE _ _ _ _ _ 

LLLL  1=0. 

KKKl=KKK  +  1 

UC  2  0  2  JJ  J  =  KKK1  ,i\l  t  I _ _ 

zzzzi  =  zzzzi*Aes(Tvs(jjjn 

IF(  ZZZZl/ AdSl  TVS(  KKK)  )  .  GE  .  1  )  I  M)=  1Nl>  +  1 

OC  203  J  J  J  =  1  »  M  1  I _ 

IF  (  A13SI  R2  (  J  JJ  )  J  .1  T.  EPS)  GC  TG  203 
KKK= JJJ 

GC  TC  204 _ 

CCMINUE 
ZZZZ1=0. 

K  K  K  1  =  KKK  +  1 _ 

GC  205  JJ J=KKKltNI  I  I  . 

ZZZZl  =  ZZZZl  +  Abi(H2 ( JJJ  I) 

1F(  ZZZZl/ARSI f<  2  ( K  KK  ) 1 .  G  E  . 1 ) I  NC  =  I  NO  +  1 
OC  2 06  JJJ=1»M  II 

I F  <  AtsSITWZl  JJJM.LT.EPS)  GC  TC  20o 

KKK= JJJ _ _ _ 

GC  TC  20  7  .  . . . 

206  CCMINUE 

2  07  ZZ7Z1=Q. _ _ 

KKKl=KKK  +  1 

DO  203  JJJ  =  KKK1 »MI  1  I 

203  ZZZZ1=ZZZZ1 ♦ A ES (TWZ (JJJ  >  ) _ 

IF (ZZZZl/ABS(TWZ( KKK) ).GE.l ) 1  ND=  I  NO  ♦  1 
i\  111  =  Mil  +  1 

_ ZZZe(L)  =  RC(MII) _ 

□  C  209  J  J  J=2  ,IMI  i  l 

209  ZZZtll  JJJ»=RC(M  II-JJJ  +  Lt  *•  T* ZZ Z d< J J J- 1 ) 

Y  N  E  «  1  1  )  =  Z  Z  Z  E  (  N  IJ  U _ 

Lit  oil) =a  l ( M  i  i ) 

OC  2  10JJJ=2  fiNI  I  I 

Z  Z  Z  ti  (  J  J  J  )  =  p  1  I  M  I  i  -  J  J  J  +  1)  »  T»  ZZZQ  <J_J  J-J  ) 
YNt*<  2)=ZZZU<  NI~1I  ) 

zzze<  i  )  =  K  2 ( M  II  ) 

OC  2  1  1  J  J  J  =  2  «M  1  1 _ 

Z  Z  ZF.  (  J  JJ  )  =K2  t  M  fi- J  jj+'u  «-  («ZZZu< JJJ-i) 
YNEW(3)=ZZZo(MII  ) 

Z  7  Z  E  1  1  ) =X VQ ( NI  l  I ) _ _ 

OC  2  1 2  JJJ  =  2  iMTV 

ZZZd(JJJ)  =  X  YO  (  I\  1 1  I  •  J  J  J  ♦  l  )  +  T*ZZZti(JJJ-l> 

ynl»I4)=z  z  z  p  mj  L 1 _ _ _ _ _ 

lLlH  (  1  )  =  XY1  (M  1  I) 

OC  21  3JJJ=2tM  I  I 

ZZZH(JJJ)=XY1(M1  1-.)JJ»1)  ♦  I-ZZZtH  JJJ-l) 

Y  IS  E  V*  <  5  )  =  Z  Z  Z  E  (  M  m 


200 

201 

202 

20  3 
204 

?  05 


k 


ZZZBI  n  =XN2  (M  I  I) 

DC  2 14 J J J  =2  »  M  1 1 

214  ZZZBI  JJJ)=XP21MI  I-JJJ  +  1)  ♦  T«Z22d(  JJJ-1) 
YNEM&)=ZZ2B(MII  ) 

Z Z Z d (  1) =  E0( M  1  I  ) 

_ DC  215JJJ-2tM  I  I _ _ 

215  ZZZB( JJJ) =EG( M  ll-JJJ+1)  +  T* ZZ ZB ( J J J- L  ) 

Y  N  E  W  I  7>=ZZZEIM  II  ) 

_ ZZZBI  l )  =E  1  (  M  I  I .) _ 

0 C  2 ldJ J J=2  *M  I  I 

216  ZZZEI  JJJ)=E1IM  Il-JJJ+l)  +  T*ZZZB( JJJ-1  ) 

_ YNEw(fl)=Z2ZI‘<MII) _ 

ZZZBI  1  )  =62  {  M  I  I  I 
DC  21 7J J J  =  2  »  N I  l  I 

217  ZZZBI  JJJ)=E2(M  II-JJJ  +  1)  »  7+ZZZ  0  (  JJJ-1  ) _ 

YKEW(91=ZZZE(M  II  > 

ZZZBI 1 ) =  T  N  P  1  ( M  II  ) 

_ DC  2 1 8 JJJ  =2  » N  I  I  I _ 

2 1  d  ZZZBI  JJJ)  =  Tf'Pl(M  II-JJJ  +  11  +  T*ZZZB(  J  J  J-  1 » 
YNEhl  10)  =  ZZZB(MI  I) 

_ ZZZBI  1)=TNP2(M  II  ) _ 

IJC  2  19 J J J=2  tNI  II 

219  ZZZBI  JJJ)  =  T.vp2(M  ll-JJJ  +  1)  +  T*  ZZ  ZB  I J  J J- 1 ) 

_ YNEM  11  )  =  ZZZB(M1  I  ) _ 

ZZZBI  L)*XK1(MI  1) 

DC  220JJJ=2  t  M  1 1  I 

220  ZZZBI  JJJ)  =XK1  I  MI  I-JJJ+1)  +  T»  Z  Z  ZB  1  J  J  J- 1 1 
YNExI  12)=ZZZB(MI  I) 

ZZZBI  L  )  s  XK2  ( M  1 1 ) 

_ DC  22  1 J J J=2  iM  1  1 _ _ 

221  ZZZBI  JJJ)=XK2(MI  I-JJJ+1)  +  T*ZZ/.B(  J  J  J- l ) 
YNEM  13  )  -  ZZZB  ( N  1 1  I  I 

_ PETLKN _ 

END 


References 


1.  J.  E.  Fromm,  "The  Time-dependent  Flow  of  an  Imcompressible  Fluid," 

Methods  in  Comput.  Physics,  _3,  New  York,  p.  346  (1964), 

2.  R.  D.  Richtmeyer  and  K.  W.  Morton,  Difference  Methods  For  Initial 
Value  Problems,  Interscience,  New  York  (1967). 

3.  0.  C.  Zienkiewicz,  The  Finite  Element  Method,  Third  Edition, 

McGraw-Hill,  London  (1977). 

4.  R.  L.  Brown,  "Investigation  of  the  Computational  Aspects  of  the 
Numerical  Solution  of  Flow  on  a  Cone."  Report  22,  Vol.  II, 

Research  Reports  of  1978  USAF-ASEE  Faculty  Research  Program, 

(1978). 

5.  S.  C.  Lubard,  and  W.  S.  Helliwell, "Calculation  of  the  Flow  on  a 
Cone  at  High  Angle  of  Attack,"  pp.  5-12,  RDA-TR-150,  R&D  Associates, 

Santa  Monica,  CA  (1973). 

6.  J.  Douglas,  and  G.  E.  Gunn, "A  General  Formulation  of  Alternating 
Direction  Methods,"  Num.  Math.,  6^,  pp.  228-253  (1964). 

7.  R.  M.  Beam,  and  R.  F.  Warming,  "An  Implicit  Finite-Difference 
Algorithm  for  Hyperbolic  Systems  in  Conservation  Law  Form," 

Journal  Computational  Physics ,  22,  pp.  87-110  (1976). 

8.  D.  W.  Peacemen,  and  H.  H.  Rachford,  "The  Numerical  Solution  of 
Parabolic  and  Elliptical  Differential  Equations,"  SIAM  Journal, 

3,  pp.  28-41  (1955). 

9.  W.  S.  Helliwell,  and  S.  C.  Lubard,  "An  Implicit  Method  for  Three- 
Dimensional  Viscous  Flow  with  Application  to  Cone  at  Angle  of 
Attack,"  Report  TR-0074(4450-64)-l,  The  Aerospace  Corporation, 

Santa  Monica,  CA  (1973). 

10.  J.  M.  Hyman,  "A  Method  of  Lines  Approach  to  the  Numerical  Solution 
of  Conservation  Laws,"  Advances  In  Computer  Methods  for  Partial 
Differential  Equations  III,  pp.  313-321,  R.  Vichnevetsky  and 

R.  S.  Stepleman  (ed.)  IMACS  (1979). 

11.  G.  Dahlquist,  "A  Special  Stability  Problem  for  Linear  Multistep 
Methods,"  BIT,  PP-  27-43  (1963). 

12.  S.  C.  Lubard  and  W.  S.  Helliwell,  "Calculation  of  the  Flow  on  a 
Cone  at  High  Angle  of  Attack,"  pp.  26-30,  RDA-TR-150,  R&D  Associates 
Santa  Monica,  CA  (1973) . 

13.  R.  M.  Beam  and  R.  F.  Warming,  "On  the  Construction  and  Application 
of  Implicit  Factored  Schemes  for  Conservation  Laws,"  SIAM-AMS 
Proceedings  of  Symposium  on  Computational  Fluid  Dynamics,  Vol.  II  (1977). 


14.  R.  L.  Brown,  "Stability  of  Sequences  Generated  by  Nonlinear 
Differential  Systems,"  Math.  Comp.,  33 ,  pp.  637-645  (1979). 

15.  G.  Dahlquist,  "G-stability  is  Equivalent  to  A-stability , "  BIT, 

18,  pp.  384-401  (1979). 

16.  N.  Rouche,  P.  Habets,  and  M.  Laloy,  Stability  Theory  by  Liapunov's 
Direct  Method,  Springer- Verlag,  New  York  (1977). 

17.  K.  R.  Kovach,  "A  Precompiler  for  Deriving  the  Time  Series  Solution 
to  Systems  of  Differential  Equations,"  M.S.  Thesis,  University 

of  Virginia,  Charlottesville,  VA  (1980). 

18.  D.  H.  Norrie,  G.  de  Vries,  The  Finite  Element  Method,  Academic 
Press,  New  York  (1973). 

19.  A.  R.  Mitchell,  R.  Wait,  The  Finite  Element  Method  in  Partial 
Differential  Equations,  John  Wiley  &  Sons,  New  York  (1977). 

20.  A.  C.  Hindmarsh,  "The  LLL  Family  of  Ordinary  Differential 
Equation  Solvers,"  UCRL-78129  (1976). 

21.  G.  D.  Byrne,  A.  C.  Hindmarsh,  "A  Polyalgorithm  for  the  Numerical 

Solution  of  Ordinary  Differential  Equations,"  ACM  TOMS ,  (1975). 


i 


DISTRIBUTION  LIST 


Copy  No. 
1-16 


17  -  18 

19 

20 


21  -  22 


Air  Force  Office  of  Scientific  Research 
Bolling  Air  Force  Base 
Washington,  D.  C.  20332 

R.  L.  Brown 

J.  M.  Ortega 

I.  A.  Fischer 

Office  of  Sponsored  Programs 

E.  H.  Pancake 
Clark  Hall 


23  RLES  Files 


0692: jt 


m 


m 


UNIVERSITY  OF  VIRGINIA 
School  of  Engineering  and  Applied  Science 

The  University  of  Virginia's  School  of  Engineering  and  Applied  Science  has  an  undergraduate 
enrollment  of  approximately  1,000  students  with  a  graduate  enrollment  of  350.  There  are  approximately 
120  faculty  members,  a  majority  of  whom  conduct  research  in  addition  to  teaching. 

Research  is  an  integral  part  of  the  educational  program  and  interests  parallel  academic  specialties. 
These  range  from  the  classical  engineering  departments  of  Chemical,  Civil,  Electrical,  and  Mechanical  to 
departments  of  Biomedical  Engineering,  Engineering  Science  and  Systems,  Materials  Science,  Nuclear 
Engineering,  and  Applied  Mathematics  and  Computer  Science.  In  addition  to  these  departments,  there  are 
interdepartmental  groups  in  the  areas  of  Automatic  Controls  and  Applied  Mechanics.  All  departments  offer 
the  doctorate;  the  Biomedical  and  Materials  Science  Departments  grant  only  graduate  degrees. 

The  School  of  Engineering  and  Applied  Science  is  an  integral  part  of  the  University  (approximately 
1,400  full-time  faculty  with  a  total  enrollment  of  about  14,000  full-time  students),  which  also  has 
professional  schools  of  Architecture,  Law,  Medicine,  Commerce,  and  Business  Administration.  In  addition, 
the  College  of  Arts  and  Sciences  houses  departments  of  Mathematics,  Physics,  Chemistry  and  others 
relevant  to  the  engineering  research  program.  This  University  community  provides  opportunities  for 
interdisciplinary  work  in  pursuit  of  the  basic  goals  of  education,  research,  and  public  service. 


I 

I 


^tkgrnmm 


